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I.  Introduction 


Most  of  the  research  done  during  the  contract  period  is  reported  in  the  “Interim  Report”  of 
March,  2002.  This  is  included  in  section  II.  Other  results,  including  a  breakthrough  that  has 
lead  to  a  new,  simpler,  fully  conservative  Vorticity  Confinement  formulation  will  be  discussed 
in  section  III.  Results  of  the  original  formulation  for  blunt  body  flow,  including  comparisons 
with  experiment,  are  described  in  section  IV.  A  new,  more  effective  way  to  treat  the  surface 
boundary  layer  that  can  be  compared  with  other  Vorticity  Confinement  separation  results, 
such  as  dynamic  stall,  is  presented  in  section  V.  A  few  of  the  publications  produced  during 
the  contract  period,  which  contain  relevant  material  refered  to  in  this  report,  are  included 
a s  appendices. 


II.  Interim  Report,  March  25,  2002 

II.  1  Introduction  to  Interim  Report 

Thin,  fixed  or  deformable  lifting  surfaces  are  very  important  components  of  many  Army 
aerodynamic  vehicles  and  devices.  These  include  rotorcraft  blades,  tail  fins,  lifting  devices 
attached  to  them  such  as  leading  edge  and  trailing  edge  flaps  and  landing  devices  such  as 
parachutes,  to  name  a  few.  It  is  well  known  that  the  development  and  incorporation  of 
any  aerodynamic  technology  such  as  described  above  requires  a  large  number  of  accurate 
simulations  to  be  performed.  Due  to  the  high  Reynolds  number,  complex  geometry  and 
aeroelastic  effects,  wind  tunnel  testing  is  very  expensive,  for  these  cases  limited  to  only 
a  small  number  of  tests,  and  in  some  cases,  not  feasible  in  current  facilities  at  full  scale. 
Hence,  it  is  necessary  that  computational  simulations  be  done.  These  must  be  efficient, 
requiring  short  set-up  and  computing  times  if  they  are  to  be  effective  in  the  engineering 
design  process.  Unfortunately,  conventional  computational  methods  fall  far  short  of  this 
requirement  for  these  flows:  Conventional  grid-based  Euler  or  viscous  (  Navier  Stokes  ) 
turbulent  modeling  methods  require  surface  fitted  or  adaptive  grids  which  often  must  be 
very  fine  near  the  surface,  and  require  far  too  much  computing  and  set-up  time  to  be  useful 
for  routine  engineering/design  applications.  This  problem  is  even  more  serious  when  the 
surface  is  deforming  so  that  grid  must  be  continually  regenerated.  Further,  conventional  grid 
based  methods  also  cannot  feasibly  be  used  when  convecting  vortical  effects  are  important 
because  of  numerical  dissipation.  Finally,  integral  methods  such  as  panel  methods,  while 
being  efficient,  cannot  treat  general  separating  cases  with  associated  vortical  effects. 

Over  the  past  several  years  we  have  developed  a  method  that  alleviates  all  of  the  above 
problems,  but  for  blunt  body  flows  .  This  method  allows  the  rapid,  simple  treatment  of 
complex  geometries,  including  modeling  of  the  relevant  effects  of  turbulent  boundary  layers, 
and  general  separation,  as  well  as  the  computation  of  convecting,  thin  vortical  regions,  with 
no  numerical  diffusion.  The  method  is  based  on  our  Vorticity  Confinement  computational 
technique.  It  utilizes  uniform  Cartesian  grids  with  “immersed  boundaries”  for  the  blunt 
body  surfaces  as  well  as  for  the  convecting  vortices.  The  basic  technique  is  described  in 
Refs.  (1-6). 
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The  reason  that  this  original  method  is  not  currently  used  for  thin  surfaces  is  that  it 
requires  that  there  be  “ghost”  points  in  the  body  to  extrapolate  the  external  flow  variables. 
This  poses  no  problem  for  the  intended  (large  class  of)  blunt  body  flows  where  at  least 
several  grid  points  can  be  included  across  the  body,  i.e.,  the  body  must  be  at  least  several 
grid  cells  thick.  This  method  has  been  validated  for  a  number  of  cases  (1_6)  and  shown  to 
be  accurate  for  these  flows. 

Currently,  an  extension  of  the  method  is  being  implemented  which  results  in  more  ac¬ 
curate  inviscid  computations  for  attached  flow  regions  (second  order  instead  of  first)  and 
utilizes  a  separate  boundary  layer  integral  computation  on  the  surfaces.  This  new  extension 
results  in  simpler  turbulent  boundary  layer  modeling  than  the  original  Confinement  method, 
since  it  utilizes  a  lower  dimensional  (surface)  grid.  For  example,  in  3-D  it  can  use  2-D  surface 
points,  which  can  be  dense,  rather  than  the  coarser  3-D  flow  grid  points  near  the  surface. 
This  new  technique  is  described  below  in  Sec.  II.2.  The  current  version  also  uses  points  in¬ 
side  the  body  and  hence  has  the  same  limitation  to  thick,  blunt  bodies  as  the  basic,  original 
Vorticity  Confinement  method.  However,  since  it  involves  solving  boundary  layer  equations 
directly  on  the  surface  and  coupling  them  to  an  outer  inviscid  solution  on  a  fixed,  uniform 
grid,  it  can  be  extended  to  treat  thin,  embedded  surfaces. 

II. 2  Development  of  Simple  Boundary  Layer  Model  Including  Sep¬ 
aration  Using  Surface  Coordinates 

During  the  performance  of  our  current  ARO  contract  involving  blunt  bodies:  “Computa¬ 
tion  of  Separating  High  Reynolds  Number  Incompressible  Flows  Using  Uniform  Cartesian 
Grids”,  a  new  extension  was  discovered  that  leads  to  far  more  efficient  computation  of  sep¬ 
arating  flows,  with  far  greater  control  over  the  boundary  layer  dynamics  than  our  original 
Vorticity  Confinement  based  method.  This  method  is  termed  “Surface  Boundary  Layer 
Model”  (SBLM).  Further,  in  addition  to  Cartesian  grids,  the  new  SBLM  method  can  be 
implemented  in  existing,  conventional  inviscid  codes  which  utilize  body  conforming  grids, 
giving  them  the  ability  to  model  boundary  layer  effects  leading  to  blunt  body  separation 
from  smooth  surfaces,  with  essentially  no  increase  in  computational  requirements. 

II.2.1  Rationale  for  SBLM 

The  aerodynamic  forces  on  surfaces  of  blunt  bodies  are  determined,  to  a  large  extent,  by 
the  zone  of  separation.  This,  in  turn,  depends  on  the  dynamics  of  the  boundary  layer  (BL), 
which  is  turbulent  for  most  flows  of  interest.  The  detailed  dynamics  can  be  solved,  without 
averaging,  by  direct  Navier  Stokes  simulation  (“DNS”),  but  only  for  very  small  regions  of 
the  BL  at  realistic  Reynolds  numbers.  Thus,  for  general  flows,  the  detailed  dynamics  cannot 
be  feasibly  computed  and  must  be  modeled  in  terms  of  some  averaged  variables.  Our  general 
approach  is  to  treat  the  spatial  resolution  of  the  BL  in  a  way  that  is  consistent  with  the 
temporal  resolution:  Since  the  BL  flow  variables  are  treated  as  averaged  over  short  periods 
of  time,  as  in  typical  “Reynolds  Averaged  Navier  Stokes”  (RANS)  treatments,  we  maintain 
that  a  simple,  consistent  model  should  involve  flow  variables  that  are  also  averaged  over 
small  regions  of  space.  Since  the  time  averaged  flow  (when  attached)  is  smooth  along  the 
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surface,  the  averaging  is  most  important  normal  to  the  surface  -  through  the  BL.  This 
amounts  to  using  variables  that  represent  BL  quantities  integrated  normally  through  the 
BL,  such  as  a  (small  number  of)  moments.  This  is  very  different  from  the  common  RANS 
approach  of  trying  to  solve  for  the  spatial  details  (in  the  normal  direction)  without  spatial 
averaging  of  some  time-averaged  profile  within  the  BL,  which  requires  very  fine  grids  which 
are  body  conforming.  Our  treatment  is  consistent  with  the  dominance  of  large  scale  coherent 
structures  in  the  BL  for  which  the  temporal  fluctuations  are  comparable  to  the  spatial,  which 
are,  in  turn,  comparable  to  the  BL  thickness^ . 

One  method  for  constructing  such  a  crude,  averaged  BL  treatment  with  a  minimal  number 
of  variables  is  to  use  the  original  Vorticity  Confinement  method(1-6).  This  results  in  thin 
vortical  BL  regions  spread  over  ~  2  grid  cells  normal  to  the  surface,  which  essentially  have 
only  a  few  degrees  of  freedom. 

II.2.2  Surface  Boundary  Layer  Model  (SBLM) 

The  subject  of  this  proposal  is  the  development  of  a  BL  model  that  can  treat  thin  lifting 
surfaces  less  than  one  grid  cell  thick.  Since  our  original  BL  scheme  requires  grid  points  inside 
the  boundary  (on  which  velocity  is  set  to  zero),  this  was  not  possible.  The  new  method, 
however,  is  far  simpler  than  the  original  Vorticity  Confinement-based  scheme,  and  much 
more  amenable  to  analysis  and  as  stated,  can  be  extended  to  thin  surfaces  extension.  In 
addition,  it  does  not  have  the  problems  of  loss  of  resolution  due  to  the  ~  2  cell  BL  thickness 
of  the  original  method  which  would  increase  the  thickness  of  a  thin  surface  by  up  to  4  cells. 
As  a  result,  since  it  has  no  displacement  thickness  errors,  the  new  method  can  even  treat 
small  features  that  have  essentially  zero  thickness.  In  addition,  the  method  can  still  be  used 
with  thick,  blunt  bodies. 

The  main  idea  is  to  combine  a  new,  simple  inviscid  immersed  boundary  treatment  with 
a  BL  model  that  “lives”  on  the  surface.  The  inviscid  solution,  which  is  smooth,  is  developed 
on  the  uniform  Cartesian  grid  outside  the  surface.  The  pressure  from  this  solution  is  extrap¬ 
olated  onto  the  separate  surface  grid  nodes  using  a  special  inviscid  flow  model.  Then,  a  lower 
dimensional  partial  differential  equation  is  solved  on  the  surface  for  the  (model)  tangential 
velocities  which  represent  the  integral  in  the  normal  direction  of  the  time  averaged  BL  veloc¬ 
ities.  The  tangential  gradient  of  the  computed  surface  tangential  velocities  is  then  used  to 
compute  a  normal  velocity  at  the  outer  edge  of  the  boundary  layer,  through  the  continuity 
equation.  In  an  adverse  pressure  gradient,  the  divergence  of  the  tangential  gradient  of  this 
tangential  velocity  leads  to  a  sudden  increase  in  normal  BL  velocity  (BL  “eruption”),  which 
initiates  separation. 

Although  it  may  seem  that  this  new  method  requires  additional  data  (a  surface  grid) 
compared  to  the  original  Vorticity  Confinement-based  BL  method,  which  only  utilized  the 
uniform  Cartesian  grid,  this  is  not  true:  the  original  surface  description  in  all  of  our  appli¬ 
cations  involves  a  set  of  surface  nodes  (usually  surface  triangles),  which  are  used  to  compute 
a  level  set  “distance”  function  on  each  Cartesian  grid  node.  This  function  then  defines  the 
surface.  Thus,  this  surface  grid  structure  is  already  required;  we  are  now  going  to  have  a 
second  use  for  it  in  the  new  BL  method. 

The  new  BL  method  is  explicitly  consistent  with  the  scenario  that  the  physical  BL  is 
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very  thin  before  separation  and  does  not  affect  the  outer,  inviscid  flow,  as  it  evolves  along 
the  surface,  driven  by  the  pressure  gradient  of  the  outer  flow  and  surface  friction,  and  that, 
when  it  “erupts”,  it  initiates  separation  in  the  outer  flow. 

First,  the  inviscid  solution  method  will  be  described,  then  the  surface  BL  model.  It  will 
be  seen  that  not  only  is  the  BL  modeled,  but,  for  the  inviscid  solution,  the  near  (sub-grid 
scale)  region  just  exterior  to  the  surface  is  also  modeled. 

II.2.2.1  Near-Surface  Inviscid  Solution  Method 

As  in  our  original  method,  we  still  discretize  and  solve  the  continuity  and  momentum 
equations  with  Vorticity  Confinement  terms  on  a  uniform  Cartesian  grid*' 1  6^.  The  main 
difference  is  that  the  velocity  on  the  grid  nodes  near  the  surface  (within  one  cell)  are  deter¬ 
mined  differently:  The  tangential  velocity  is  determined  by  Surface  Vorticity  Confinement, 
which  effectively  enforces  zero  vorticity  in  the  grid  cells  exterior  to  the  surface  region.  The 
inviscid  near-surface  “model”  referred  to  above  then  specifies  the  normal  velocity  at  nodes 
in  the  near-region  (see  Fig.  1). 

Even  though  the  BL  does  not  result  from  a  balance  of  diffusion  and  confinement,  as 
it  did  in  the  original  Confinement  based  method,  Vorticity  Confinement  is  still  vital  here: 
Numerical  errors  normally  create  vorticity  which,  without  Confinement,  will  diffuse  and 
convect  away  from  the  boundary  region,  contaminating  the  outer  flow  field  with  large  errors: 
Confinement  prevents  this  by  convecting  erroneous  vorticity  back  into  the  boundary  surface, 
and  results  in  a  well-defined  outer  inviscid  flow. 


Figure  1:  Computational  Method 
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The  near-surface  inviscid  model  specifies  that  Vsqs  is  independent  of  normal  distance, 
n,  in  the  near  boundary  region,  where  Vs  is  the  derivative  along  the  surface  (in  our  initial 
work,  we  approximate  the  surface  as  a  flat  plate).  This  inviscid  model  is  reasonable  because 
the  outer  vorticity  is  zero  before  separation,  implying 

&nQs  ^ sQn  0 

where  the  normal  velocity,  qn,  is  zero  on  the  boundary,  and  its  n  derivative  varies  smoothly 
along  the  surface.  Thus  qn  is  0(h)  (grid  cell  size)  except  at  separation.  Accordingly,  dnqs  is 
0(h)  and  qs  only  varies  by  0(h)  in  the  near  surface  “inner”  modeled  region,  of  thickness  h. 
Then,  Vs  •  qs  can  be  taken  to  be  independent  of  n  to  0(h).  At  a  point  n  above  the  surface 
(n  ~  0(h)),  then,  we  can  model  qn  by  integrating 

dnqn  =  -V.  •  qs 


This  gives 

qn  ~  —n(V s  •  <fs)  |n=0  +  <?n 

where  qbn  represents  the  effect  of  the  boundary  layer  model,  and  is  described  in  Sec.  II.2.2.2. 
The  value  of  qn  is  enforced  at  nodes  just  outside  the  surface  (see  Fig.  1). 

Results  of  an  inviscid  computation  (without,  of  course,  separation)  of  a  vortex  in  a 
straight  channel  (not  aligned  with  the  grid)  using  the  new  method  are  presented  in  Fig.  2. 


Figure  2:  Velocity  vector  of  an  inviscid  computation  of  a  vortex  in  a  straight  channel 

There  are  many  opportunities  to  include  effects  such  as  surface  curvature  here,  as  well  as 
other  near  surface  inviscid  models.  For  example,  for  long,  thin  lifting  surfaces,  a  lifting-line 
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model  for  main  and  tail  rotor  blades  with  sub-grid  scale  chord  was  recently  implemented 
as  a  “near  surface  inviscid  model”  for  a  project  for  United  Technology  Research  Corp.^. 
There,  only  a  uniform  Cartesian  grid  was  used  which  the  lifting  lines  moved  through.  This 
project  resulted  in  a  code  able  to  predict  flow  over  an  entire  rotorcraft  including  fuselage, 
empennage  including  tail  and,  as  described,  main  and  tail  rotors.  Here,  the  rotor  surface 
was  a  blade  with  chord  assumed  so  small  that  it  was  less  than  a  grid  cell  and,  as  a  result,  the 
blade  could  be  approximated  by  a  vortex  filament  with  specified  motion.  (This  project  really 
was  a  preliminary,  simplified  application  of  the  inviscid  part  of  the  model  to  thin  surfaces). 

II.2.2.2  Boundary  Layer  Model 

The  new  BL  model  is  effected  directly  on  nodes  on  the  surface.  This  model  is  easy  to 
analyze  and  modify  with  additional  terms  since  it  exists  on  a  lower  dimensional  surface  and 
is  driven  by  a  pressure  field  extrapolated  from  a  smooth,  exterior  inviscid  solution.  The 
equations  to  be  solved  along  the  surface  are: 

dt<ts  =  -9$  •  V<?1  -  Vs(pe/p)  +  fs  (1) 

where  (fs  is  the  model  tangential  velocity  for  the  boundary  layer:  it  is  single  valued  at 
each  point  and  represents  a  model  for  the  velocity  averaged  through  the  layer  in  the  normal 
direction  as  well  as  averaged  in  time.  Also,  Vs(pe/p)  is  the  tangential  derivative  of  the 
inviscid,  near  surface  outer  pressure  field  extrapolated  to  the  surface,  divided  by  the  density, 
p.  The  surface  effect,  fs,  can  be  represented  by  a  number  of  forms.  This  is  essentially 
our  “turbulent  boundary  layer  model”.  Without  it,  there  is  no  separation  and  the  surface 
velocity  equals  the  near-surface  inviscid  velocity.  A  simple,  single  parameter  “drag”  model 
that  we  are  using  initially  is 

fs  =  -Ml 

This  is  one  of  the  simplest  possible  forms  and  is  an  obvious  first  choice.  This  model  was 
studied  for  flow  over  a  2-D  cylinder  induced  by  a  pair  of  convecting  vortices  as  they  ap¬ 
proached  a  cylinder^.  There,  it  was  shown  to  lead  to  a  singularity  in  Vs^.  This  is  expected 
since  eqn.  (1)  is  just  Burger’s  equation  with  an  imposed  forcing  function  and  it  is  known^ 


for  tc,  c  constants.  That  is,  the  gradient  diverges  at  a  particular  position  and  time. 

The  last  part  of  the  BL  model  simulates  the  effect  of  incompressibility.  It  represents,  in  a 
simple  way,  what  happens  in  more  detailed  computations  of  a  separating  boundary  layer(10); 
when  the  tangential  gradient  of  the  tangential  velocity  diverges,  the  gradient  of  the  normal 
velocity  diverges,  since  incompressibility  implies 

a„^  =  -vs-q1 

Taking  $  to  be  a  mean  velocity  through  the  BL  with  assumed  thickness  S  <  h,  with  qn  =  0 
at  the  surface,  this  results  in 

s 

9n  =  -/v-9ldn 
o 
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qbn  =  -Ws-Q ts 

This  normal  velocity  is  added  to  the  model  normal  velocity  on  the  outer  grid  nodes  for  the 
inviscid  near-boundary  solution  described  above  in  Sec.  II.2.2.1. 

During  most  of  the  time  qbn  will  be  small  and  have  little  effect,  since  V  •  q%  is  0(1)  and 
S  is  small.  However,  when  the  gradient  diverges  qbn  will  become  large  and  result  in  a  large 
qbn  and,  hence,  a  BL  eruption.  While  the  detailed  BL  computation  of  Ref.  (10)  had  to  be 
stopped  at  this  point  because  it  did  not  treat  the  outer  flow,  ours  can  continue.  Preliminary 
results  of  such  an  eruption  (prescribed  qbn)  have  shown  that  it  leads  to  separating  vortices(6). 

II.2.3  Current  Work 

During  this  final  phase  of  our  current  ARO  project,  we  are  solving  vortex-induced  separation 
in  2-D,  and  3-D  flows  over  circular  cylinder,  ellipsoid  and  sphere,  with  different  grid  sizes, 
in  order  to  calibrate  the  BL  model  and  determine  the  best  values  for  the  coefficients  ft  and 
S  and  whether  other  terms  (such  as  tangential  derivatives  or  surface  curvature  dependence) 
are  necessary,  or  whether  //  and  S  should,  themselves,  be  solutions  to  a  surface-transport 
equation. 

Fig.  3  show  the  comparison  of  flow  over  cylinder  with  exact  solution.  Turbulent  boundary 
layer  effects  can  be  modeled  as  shown  in  Fig.  4.  Fig.  5  and  Fig.  6  show  the  velocity  vector 
of  the  inviscid  flow  over  an  ellipsoid  and  a  sphere  respectively. 


(a)  Inviscid  cylinder  with  SBLM  (b)  Exact  solution 

Figure  3:  Comparison  of  Velocity  Vector  of  Flow  over  Cylinder 
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(a)  fi  =  0.01  (b)  M  =  0.005  (with  delayed  separation) 

Figure  4:  Vorticity  Contour  of  Vortex  Shedding  (Model  of  Turbulent  B.L.  Effects) 


~  _ 

.  --  - - —  - - — - — 


(a)  Side  View  (b)  Top  View 

Figure  5:  Velocity  Vector  of  Inviscid  Flow  over  Ellipsoid 
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Figure  6:  Velocity  Vector  of  Inviscid  Flow  over  Sphere 


II.3  Development  of  Surface  Boundary  Layer  Model  for  Thin  Sur- 
faces 

As  explained  above,  computation  of  flow  over  general,  thin,  fixed  surfaces  is  currently  not 
feasible,  with  conventional  methods,  for  rapid  engineering  computations:  For  general,  sepa¬ 
rating  flow  the  very  fine  body  fitted  grids  necessary  for  conventional  viscous  computations 
cannot  be  feasibly  generated  because  the  surfaces  generally  are  thin  and,  further,  can  move 
in  response  to  aerodynamic  or  actuator  forces,  requiring  continual  re-generation.  Further, 
current  pde-based  turbulent  models  are  too  expensive  to  routinely  run  for  changing,  com¬ 
plex  geometries.  Our  SBLM  approach  totally  alleviates  these  problems  while  not  requiring 
any  more  turbulence  modeling  parameters  than  conventional  pde-based  turbulent  modeling 
schemes.  This  is  because  only  a  regular,  coarse  inviscid-size  Cartesian  grid,  is  used  together 
with  a  separate  surface  grid,  on  which  a  simplified  turbulence  model  is  solved,  so  that  the 
computation  is  very  efficient. 

As  explained  above,  a  first,  simplified  example  of  this  approach  has  already  been  demon¬ 
strated  for  United  Technology  Research  Corp/8)  using  a  flexible  lifting  line  model  for  the 
rotor  blades  (both  main  and  tail)  together  with  table  look-up,  coupled  to  the  main  flow 
field  computation.  Of  course,  the  blades  were  not  treated  as  surfaces  and  no  boundary  layer 
model  was  employed. 

A  next  step,  suggested  by  Frank  Caradonna,  involved  computing  inviscid  flow  over  a 
zero  thickness  flat  plate  at  different  angles  of  attack  in  2-D.  Basically,  this  demonstrates  the 


accuracy  of  the  coupling  of  the  near  surface  inviscid  model  to  both  the  outer  flow  field  and 
the  surface  equation  solver  (see  Fig.  1).  Here,  the  surface  equations  only  involved  an  inviscid 
pressure  computation  and  integration  (to  get  the  lift  coefficient).  Even  though  no  turbulence 
model  was  involved,  the  resolution  requirements  of  this  problem  are  very  demanding  since 
the  pressure  is  singular  at  the  leading  edge.  Our  surface  grid  approach  allows  us  to  resolve 
this  case  since  the  surface  grid  can  be  far  finer  than  the  (higher  dimensional)  regular  field 
grid.  As  in  the  above  lifting  line  approach,  this  study  was  really  a  preliminary  test  of  the 
concept  of  the  work  proposed  here. 

Besides  demonstrating  an  important  aspect  of  our  approach,  this  case  is  extremely  im¬ 
portant  for  practical  rotorcraft  computations:  even  with  fixed,  non-deforming  surfaces,  there 
are  always  important  lifting  surfaces  such  as  tails  that  cannot  be  treated  as  lifting  lines  since, 
unlike  rotor  blades,  they  have  a  small  aspect  ratio.  On  the  other  hand,  they  are  too  thin 
for  our  original  Vorticity  Confinement  treatment,  which  would  require  that  the  outer  grid 
be  fine  enough  to  have  cells  within  the  thickness  of  the  surface,  resulting  in  fine  grids  with 
very  large  numbers  of  cells.  This  was  adequate  for  our  original  blunt-body  based  surface 
treatment  but  not  for  the  present  proposed  project,  and  was  the  main  drive  for  developing 
the  new,  thin-surface  method.  Further,  the  assumption  of  zero  thickness  here  forms  a  good 
basis,  since  effects  of  finite,  small  thickness  for  thin  lifting  surfaces  are  usually  small  and,  if 
needed,  are  best  treated  as  a  perturbation  of  our  simple,  efficient  method. 

The  basic  idea  for  the  flat  plate  is  to  use  a  closed-form  solution  in  a  uniform  free  stream 
as  a  near-surface  inviscid  model.  This  is  analogous  to  a  Biot  Savart  model  for  the  lifting 
line  model  or  the  constant  tangential  velocity  and  linearly  increasing  normal  velocity  model 
for  the  infinite  flat  surface  described  in  Sec.  II.2.  Thus,  the  lifting  line  corresponds  to  zero 
chord,  the  short  plate  here  to  finite  chord  and  the  flat  surface  of  Sec.  II.2  to  infinite  chord. 

Basically,  in  all  of  the  cases,  near  the  surface  the  flow  is  assumed  to  be  the  sum  of  a  local 
uniform  free  stream  (to  be  determined  as  part  of  the  coupled  computation)  and  a  model 
term.  In  all  three  cases,  the  main  point  is  to  extrapolate  the  outer  velocity  into  the  near¬ 
surface  region  and  compute  an  effective  free  stream  velocity,  and  then  to  use  this  to  both 
get  the  effect  of  the  surface  on  the  outer  flow,  and  the  pressure  on  the  surface.  Once  this 
surface  pressure  is  computed,  our  viscous  model  which  computes  separating  flow,  described 
in  Sec.  II.2,  can  be  used. 

For  the  flat  plate,  an  exact,  analytic  solution^11)  was  used  as  the  starting  point.  This 
is  shown  in  Fig.  7,  where  V  is  the  (local)  free  stream  velocity  magnitude  and  a  the  angle 
between  local  free  stream  and  the  flat  plate.  The  equation  shown  gives  the  normal  ( q„ )  and 
tangential  (qs)  velocities  at  a  nearby  point,  P.  During  the  overall  computation,  these  values 
are  computed  at  (outer,  Cartesian)  grid  points  near  the  flat  plate  for  cell  nodes  where  the  cell 
intersects  the  plate.  Conventional  difference  equations,  together  with  Vorticity  Confinement 
are  used  elsewhere.  Matching  the  solutions  determines  the  local  V  and  9.  This  is  then  used 
in  the  equation  of  Fig.  7  to  determine  the  velocity  at  points  on  the  surface.  The  resulting 
pressure  is  then  computed  and  integrated  to  get  C).  For  the  computations,  a  256  x  256, 
uniform  Cartesian  outer  grid  was  used  and  a  1000  uniform  point  surface  grid.  Computed 
Cp  values  on  the  surface  for  angles  of  attack  (a)  of  1°,  5°  and  10°  are  presented  in  Fig.  8. 
They  are  exact  since  the  closed  form  solution  was  used  and  the  iterated  coupling  technique 
determined  the  correct  local  free-stream  flow  (the  local  flow  near  the  surface,  here  denoted 
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“free  stream”  is,  of  course,  not  input  to  the  model  but  must  be  determined  since,  in  general, 
it  is  not  known  but  is  part  of  the  flow  computation).  As  a  final  test,  the  computed  surface 
pressure  was  numerically  integrated  to  get  Ci  as  a  function  of  a.  This  is  plotted  in  Fig.  9.  It 
can  be  seen  that  it  is  very  close  to  the  theoretical  2na  line.  This  agreement  is  very  significant 
and  represents  a  breakthrough  because  of  the  singular  nature  of  the  pressure.  It  should  be 
mentioned  that  this  solution  is  only  meant  to  be  a  demonstration  of  the  matching  technique. 
The  use  of  this  analytic  closed  form  local  2-D  solution,  as  above,  should  be  useful  for  long 
thin  rotor  blades,  including  surface  effects:  It  will  be  a  great  improvement  over  lifting  line 
schemes,  while  being  just  as  fast.  The  actual  analytic  form  will  not  be  used  for  short  aspect 
ratio  surfaces  such  as  fins.  These  will  be  treated  by  the  flat-surface  technique  described 
above  (for  channel  flow). 


O=(@-0)/2 

R=  (R7R)“ 

qs =FZ?smasin0+F'cosa 
qn=VRsmacas9 


Figure  7:  Velocity  Formula  for  Flat  Plate 


III.  New  Vorticity  Confinement  Formulation 

In  spring,  2002,  we  discovered  a  radically  different  Vorticity  Confinement  formulation.  This 
new  formulation  is  completely  conservative.  The  original  formulation  had  (small)  effects 
due  to  non-conservation:  There  could  be  small  deviations  in  the  trajectory  of  vortices. 
These  could  be  corrected  by  adding  terms  to  the  original  formulation  to  make  it  fully 
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Figure  8:  Cp  Distribution  at  Different  Angles  of  Attack 


Figure  9:  Ci  Versus  Angle  of  Attack 
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conservative^12);  but  the  result  was  somewhat  complicated.  The  new  form  is  not  only  sim¬ 
pler,  but  is  amenable  to  analysis.  It  can  also  be  used  directly  for  “scalar  confinement”, 
which  is  important  for  convection  of  passive  scalars.  The  reason  for  the  simplicity  is  that 
it  only  involves  a  simple  difference  formulation  with  no  extra  logic:  A  gradient  of  vorticity 
magnitude  with  upwind  bias  was  required  in  the  original  form. 

A  detailed  rationale  is  given  in  Sec.l  of  Ref.  (13),  included  as  Appendix  I.  Results  for 
both  convecting  scalars  and  vortices  in  2-D  incompressible  flow  are  given  in  Sec. 2  of  that 
paper. 

For  concentrated,  convecting  passive  scalars,  a  new  result  is  proven  -  a  type  of  “Ehrenfest 
Theorem”  which  shows  that  the  centroid  of  the  scalar  feature  convects  with  the  mean  velocity, 
defined  as  an  average  velocity  over  the  feature  weighted  by  the  amplitude  of  the  scalar. 
Additional  discussions  are  presented  in  Sections  1.5  and  1.6  of  Ref.(14),  included  as  Appendix 

II. 

IV.  Validation  Studies  -  Original  Vorticity  Confinement 

Three  computations  have  recently  been  done  to  validate  the  original  Vorticity  Confinement 
method  for  realistic  3-D  flows. 

IV.  1 

The  first  involves  flow  over  a  rotorcraft  landing  ship.  There,  the  vortex  created  at  the 
leading  edge  corner  of  the  deck  in  a  even  moderate  cross-wind  can  convect  over  the  deck, 
remaining  strong  and  concentrated.  This  creates  hazards  for  the  landing  rotorcraft.  Results 
are  presented  in  Ref. (15)  and  Ref. (16)  (included  as  Appendix  III). 

In  section  3.2  of  the  latter,  it  can  be  seen,  first,  that  Vorticity  Confinement  was  necessary 
even  to  get  reasonable  results,  with  the  computational  grid  used.  Second,  it  can  be  seen  that 
the  vortex  position  and  induced  velocity  are  very  close  to  that  measured  experimentally. 

IV.2 

Two  canonical  3-D  turbulent  blunt  body  flows  were  computed  and  compared  with  exper¬ 
iment.  These  involved  flow  over  a  circular  cylinder  at  a  Reynolds  number  of  3,900  and 
a  square  cylinder  at  a  Reynolds  number  of  21,400.  Results  of  the  first  were  presented  in 
Ref.  (15)  and  Ref.  (14)  (included  as  Appendix  III).  The  square  cylinder  results  were  also 
presented  in  the  latter. 

The  main  goal  of  this  study  was  to  investigate  the  promise  of  Vorticity  Confinement  as  a 
way  to  effectively  and  economically  model  the  small  scale  vortices  in  a  turbulent  wake  flow. 
This  was  a  very  different  case  than  previous  trials  which  involved  isolated,  concentrated 
vortices  and  surface  shear  layers,  since  a  large  range  of  scales  was  involved.  This  study, 
then,  evaluated  Vorticity  Confinement  as  a  very  simple,  effective  alternative  Large  Eddy 
Simulation  (LES)  method.  Only  a  simple  parameter,  the  confinement  multiplier,  “e”  was 
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used  to  model  Reynolds  number  effect,  which  was  held  constant  throughout  the  field.  As  e 
was  increased,  the  small  scales  and  velocity  fluctuations  became  more  energetic. 

In  Sec.  3  of  Ref.  (16)  (Appendix  III),  the  changes  in  the  rms  fluctuating  velocity  profiles 
can  be  seen  as  e  is  varied.  For  both  circular  and  square  cylinders,  a  value  of  t  was  easily  found 
such  that  both  mean  and  rms  velocity  profiles  matched  experiment  extremely  closely.  We  are 
currently  searching  for  measurements  at  other  Reynolds  numbers  to  calibrate  e.  It  surprised 
us  that  only  the  single  confinement  term,  with  one  parameter  seemed  to  serve  as  an  effective 
LES  model.  We  believe  that  our  ability  to  capture  small  scale  vorticies  over  only  ~  2  grid 
cells  gave  us  a  large  enough  range  of  scales,  even  on  the  moderate  grid  used,  to  accurately 
capture  the  larger  scale  turbulent  dynamics.  This  is  possible  because  we  use  an  inherently 
discrete,  nonlinear  model  directly  on  the  grid.  This  is  very  different  from  conventional  LES 
schemes  which  attempt  to  model  the  small  scale  structures  with  model  partial  differential 
equation  (pde)  terms.  These  pde’s  must  then  be  resolved  by  finite  different  schemes,  which 
require  many  more  grid  points. 


V.  Extended  Surface  Boundary  Layer  Treatment 
VJ 

The  main  part  of  the  current  project  involved  developing  a  computational  method  for  blunt 
bodies  in  incompressible  flow  where  the  body  could  be  simply  represented  by  surface  coor¬ 
dinates  and  “immersed”  in  a  uniform,  coarse  Cartesian  grid. 

The  first  part  involved  the  outer,  inviscid  flow.  This  was  accomplished  in  a  simple  way 
using  a  level  set  description  of  the  surface,  and  Vorticity  Confinement  to  remove  interpolation 
errors  (which  show  up  as  numerical  vorticity)  in  the  region  near  the  surface.  The  new 
feature  here  is  that  an  inviscid  “model”  was  used  near  the  surface  eliminating  the  first 
order  numerical  boundary  layer  created  with  the  original  method,  resulting  in  a  much  more 
accurate  solution. 

The  second  part  of  the  project  involved  viscous  separation.  This  was  readily  treated  in 
an  approximate  way  by  the  original  Vorticity  Confinement  method  with  no-slip  boundary 
conditions  at  the  surface.  To  increase  the  accuracy  and  allow  additional  terms  to  be  added  to 
create  an  effective  turbulent  boundary  layer  model  that  could  be  easily  calibrated  a  “Surface 
Boundary  Layer  Model”  (SBLM)  was  developed  that  was  matched  to  the  new  inviscid  near¬ 
surface  “inviscid  model”  and  did  not  depend  on  the  orientation  of  the  surface  with  respect 
to  the  main  Cartesian  grid.  This  allowed  the  modeling  of  the  boundary  layer  separation 
in  a  grid-independent  way.  Preliminary  results  using  a  single  surface  velocity  equation  are 
reported  in  Sec.  II. 

Frank  Caradonna  (17)  brought  up  the  point  that  good  results  have  been  obtained  for 
surface  separation  using  Vorticity  Confinement  in  this  case  for  dynamic  stall^18^  -  a  difficult 
problem.  This  study  used  a  body  fitted  grid  where  the  boundary  layer  was  confined  by 
the  Vorticity  Confinement  to  ~  2  cells  near  the  surface  -  effectively  forming  an  SBLM.  Due 
to  the  body  conforming  grid  there  it  was,  of  course,  independent  of  surfce  orientation.  He 
suggested  trying  to  use  such  an  approach  (only  very  near  the  surface),  since  it  had  been 
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shown  to  be  effective. 

For  this  reason,  we  recently  have  implemented  this  approach.  Even  though  it  uses  a 
“body  fitted  grid”,  this  is  trivial  to  create  since  it  is  confined  to  the  near  surface  region  and 
only  requires  extending  the  surface  defining  coordinates  by  a  small  amount.  The  solution 
on  this  grid  is  then  matched  to  our  “outer”  inviscid  solver.  We  are  finding  that  this  method 
is  much  less  sensitive  than  our  original,  “single  surface  variable”  method,  though  it  still 
is,  effectively,  a  surface  method,  but  with  2-3D  velocity  variables  (grid  planes).  It  allows 
independent  verification/testing  of  the  surface  model,  include  post-separation.  It  also  allows 
contact  with  body-fitted  dynamic  stall  work. 

Initial  results  (using  a  much  finer  near-body  grid  than  necessary)  for  a  circular  cylinder 
are  presented  in  Fig.  10  and  11  for  the  flow  on  the  inner  planes.  Note  that  just  the  first 
2-3  planes  take  part  in  initial  separation  -  this  is  the  “Extended”  SBLM.  The  full  outer  grid 
flow  resulting  from  the  separation  is  shown  in  Fig.  12.  The  diameter  of  the  circular  cylinder 
was  16  unit  grid  cells.  The  outer  grid  has  128  x  128  unit  grid  cells  and  the  inner  grid  has 
20 (radius)  x  50 (circumferential)  grid  cells  (with  the  cell  size  dr  =  0.1  and  dd  =  27r/50). 
The  time  step  for  outer  and  inner  region  are  0.2  and  0.02  respectively.  Results  shown  was 
after  100  time  steps,  fi  for  the  outer  and  inner  region  are  0.75  and  0.05  respectively  and  e 
for  the  outer  and  inner  region  are  1.125  and  0  respectively. 
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Figure  10:  Velocity  Vector  of  the  Inner  Region  for  Flow  Over  Circular  Cylinder  (without 
interaction  with  outer  region) 
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Figure  11:  Velocity  Vector  of  the  Inner  Region  for  Flow  Over  Circular  Cylinder 


Figure  12:  Velocity  Vector  for  the  Flow  Over  Circular  Cylinder 


V.2  Conclusion  of  Section 

An  extended  SBLM  approach  was  shown  to  be  effective.  Now,  additional  model  terms  can 
be  added  in  a  robust  -  grid  independent  way  -  while  still  immersing  the  surface  in  a  simple, 
uniform  Cartesian  grid.  Additional,  more  complex  configurations,  including  a  3-D  prolate 
ellipsoid  are  being  studied. 
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Abstract 


A  new  computational  method  is  described  that  efficiently  treats  thin  features  in 
multi-dimensional  incompressible  fluid  flow,  such  as  vortices  and  streams  of  passive 
scalars,  and  convects  them  over  long  distances  with  no  spreading  due  to  numerical  errors. 
Outside  the  feature,  where  the  flow  is  irrotational  or  the  scalar  vanishes,  conventional 
discretized  finite  difference  fluid  dynamic  equations  are  solved.  The  feature  is  treated  as  a 
type  of  weak  solution  and,  within  the  feature,  a  nonlinear  difference  equation,  as  opposed 
to  finite  difference  equation,  is  solved  that  does  not  necessarily  represent  a  Taylor 
expansion  discretization  of  a  simple  partial  differential  equation  (pde).  The  approach  is 
similar  to  shock  capturing,  where  conservation  laws  are  satisfied,  so  that  integral 
quantities  such  as  total  amplitude  and  centroid  motion  are  accurately  computed  for  the 
feature.  A  more  general  approach  is  needed,  however,  than  for  conventional  shock 
capturing,  which  are  basically  1-D  phenomena.  Basically,  we  treat  the  features  as  multi¬ 
dimensional  nonlinear  discrete  solitary  waves  that  live  on  the  computational  lattice. 
These  obey  a  “confinement  ”  relation  that  is  a  rational  generalization  to  multiple 
dimensions  of  1  -D  discontinuity  capturing  schemes. 

We  have  already  been  working  for  several  years  with  a  method  -  “Vorticity 
Confinement”  that  creates  solitary  waves  on  the  lattice  over  only  2-3  grid  cells  and 
efficiently  solves  these  problems,  and  have  generated  a  large  number  of  results  for  basic 
and  complex,  realistic  flows.  The  main  objective  of  this  paper  is  to  introduce  a 
completely  new  formulation,  that,  compared  to  the  original  formulation,  is  simpler, 
allows  more  detailed  analysis,  and  results  in  more  compact  structures  for  both  vortices 
and  convecting  scalars. 

First,  some  practical  motivations  are  described.  Then,  a  short  critique  of 
conventional  methods  for  these  problems  is  given.  The  basic  new  method  is  then 
described.  Some  analysis  of  the  new  method  and  initial  results  are  finally  presented. 
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1.  Introduction 


A  new  computational  method  is  described  that  efficiently  treats  thin  features  in 
multi-dimensional  incompressible  fluid  flow,  such  as  vortices  and  streams  of  passive 
scalars,  and  convects  them  over  long  distances  with  no  spreading  due  to  numerical  errors. 
Outside  the  feature,  where  the  flow  is  irrotational  or  the  scalar  vanishes,  conventional 
discretized  finite  difference  fluid  dynamic  equations  are  solved.  The  feature  is  treated  as  a 
type  of  weak  solution  and,  within  the  feature,  a  nonlinear  difference  equation,  as  opposed 
to  finite  difference  equation,  is  solved  that  does  not  necessarily  represent  a  Taylor 
expansion  discretization  of  a  simple  partial  differential  equation  (pde).  The  approach  is 
similar  to  shock  capturing  [1],  where  conservation  laws  are  satisfied,  so  that  integral 
quantities  such  as  total  amplitude  and  centroid  motion  are  accurately  computed  for  the 
feature.  A  more  general  approach  is  needed,  however,  than  for  conventional  shock 
capturing,  which  are  basically  1-D  phenomena.  Basically,  we  treat  the  features  as  multi¬ 
dimensional  nonlinear  discrete  solitary  waves  that  live  on  the  computational  lattice. 
These  obey  a  “confinement  ”  relation  that  is  a  rational  generalization  to  multiple 
dimensions  of  1-D  discontinuity  capturing  schemes  (this  is  elaborated  on  in  the 
Appendix). 

Differences,  compared  to  conventional  shock  capturing,  are  that: 

First,  unlike  shocks,  characteristics  do  not  point  into  the  feature,  and  extra  terms 
must  be  designed  to  prevent  it  from  spreading  due  to  numerical  effects  in  the  convection. 
(Harten  [2]  developed  such  a  scheme,  but  for  contact  discontinuities  in  1-D  compressible 
flow.) 

Second,  shocks  are  basically  one  dimensional  discontinuities,  unlike  vortex 
filaments  or  thin  streams  of  passive  scalars,  which  are  intrinsically  multi-dimensional:  A 
concatenation  of  1-D  “capturing”  operators  along  separate  axes  will  not,  generally,  give 
smooth  solutions.  Due  to  the  multidimensional  nature,  it  seems  necessary  to  pay  attention 
to  the  structure  within  the  feature,  even  though  it  is  sampled  on  only  a  few  grid  cells  in 
the  cross-section. 

We  have  already  been  working  for  several  years  with  a  method  -  “Vorticity 
Confinement”  that  efficiently  creates  solitary  waves  on  the  lattice  and  efficiently  solves 
these  problems,  and  have  generated  a  large  number  of  results  for  basic  and  complex, 
realistic  flows  [3-12].  The  main  objective  of  this  paper  is  to  introduce  a  completely  new 
formulation,  that,  compared  to  the  original  formulation,  is  simpler,  allows  more  detailed 
analysis,  and  results  in  more  compact  structures  for  both  vortices  and  convecting  scalars. 

First,  some  practical  motivations  are  described.  Then,  a  short  critique  of 
conventional  methods  for  these  problems  is  given.  The  basic  new  method  is  then 
described.  Some  analysis  of  the  new  method  and  initial  results  are  finally  presented. 


1.1  Motivation 


There  are  many  important  fluid  dynamic  problems  where  thin,  concentrated  flow 
features  must  be  numerically  convected  over  long  distances.  Examples  include  vortices 
shed  from  forebodies  and  lifting  surfaces  of  aircraft  ,  rotorcraft  and  submarines,  and  thin 
streams  of  contaminants  convecting  in  ambient  flow.  Often,  for  these  cases,  the  main 
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interest  is  the  integrated  amplitude  through  the  feature  and  the  motion  of  the  centroid, 
rather  than  the  details  of  the  internal  structure.  In  general,  these  features  can  originate  in 
many  places,  reattach,  merge,  and  have  complex  topology.  Accordingly,  we  consider 
Eulerian  methods  where  very  general  topologies  can  be  treated,  as  opposed  to  Lagrangian 
“particle  tracking”  methods. 


1.2  Current  Methods 


Conventional  Eulerian  approaches  to  the  convection  problem  for  incompressible 
flow  involve  formulating  governing  pde’s,  discretizing  them  and  solving  them  as 
accurately  as  possible  on  feasible  computational  grids,  assuming  only  strong  solutions. 
For  smooth,  non-thin  features,  these  methods  are  well  known  to  converge  to  the  correct 
solution  as  the  number  of  points  across  the  feature  in  a  typical  direction,  N,  becomes 
large:  Error  estimates  are  asymptotic  in  N.  For  accurate  solutions,  even  higher  order, 
complex  discretization  methods  typically  require  N  to  be  at  least  ~8  or  10  so  that  the  error 
obeys  the  large  N  estimate.  Even  then,  it  is  well  known  that  solutions  degrade  over  long 
convection  distances  (many  time  steps).  As  a  result,  although  conventional  methods  may 
be  efficient,  for  example,  for  low  Reynolds  number  flow  problems  with  large  smooth 
vortical  regions,  they  are  inefficient  (or  not  even  feasible)  for  thin  aircraft  trailing  vortices 
convecting  over  long  distances. 

It  appears  that  adaptive,  unstructured  grids  cannot  improve  the  resolution 
significantly  for  realistic  problems  with  many  thin  features.  Also,  even  for  small  numbers 
of  time  dependent  features  these  methods  are  very  expensive,  compared  to  uniform  grid, 
structured  methods. 

For  these  reasons,  for  the  problems  considered,  it  is  important  to  have  only  very 
few  (2  or  3)  grid  points  across  each  feature  if  it  is  very  thin,  and  to  convect  it  with  no 
numerical  spreading.  This  small  number  of  grid  points  is  consistent  with  the  stated  desire 
to  only  compute  a  few  integral  quantities  across  the  feature,  such  as  total  amplitude  and 
centroid  position.  Then,  the  difference  scheme  can,  effectively,  serve  as  a  simple, 
implicit  “solitary  wave”  model  for  the  internal  structure.  However,  as  explained  above, 
this  is  not  feasible  with  conventional  pde-based  methods.  As  an  example  of  the 
limitations  of  conventional  methods,  if  we  are  limited  to,  say,  a  (128) 3  computation  and 
use  a  conventional  method  which  requires  ~8  cells  to  avoid  excessive  numerical 
spreading,  it  may  then  be  reasonable  to  say  that  the  smallest  “computational  scale”,  A,  is 
8/128  or  1/16,  rather  than  the  desired  scale  of  ~  2/128  or  1/64.  Scales  much  smaller  than 
1/16  will  be  strongly  damped  and  features  will  quickly  spread  due  to  numerical  effects  if 
a  pde  is  to  be  resolved.  This  difficulty  shows  up  as  a  limit  on  the  effective  number  of 
degrees  of  freedom  of  the  simulation,  which  is  reduced  by  a  factor  of  ~43,  if  a 
conventional  method  is  compared  to  one  with  an  actual  2  cell  small  scale  resolution. 


2.  Vorticitv  Confinement  Approach 

As  stated,  over  the  last  several  years  a  new  method  was  developed  to  treat 
problems  with  thin,  convecting  vortices  [3-12].  This  method,  termed  “Vorticity 
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Confinement”,  is  intrinsically  multidimensional  and,  although  equally  applicable  to  thin 
vortices  and  (in  modified  form),  thin  streams  of  passive  scalars,  almost  all  of  its  uses 
have  involved  vorticity.  First  (through  Sec.  2.1),  a  section  of  ref.  [3]  (with  some  changes) 
will  be  included  that  describes  the  basic  features  of  Vorticity  Confinement.  Then,  the  new 
formulation  will  be  introduced,  and  new  analyses  and  results  involving  the  application  of 
the  method  to  convecting  vortices  and  scalars  presented. 

Vorticity  Confinement  is  a  method  to  preserve  convecting,  concentrated  vorticity 
on  a  coarse  grid;  the  vortical  structures  can  represent  boundary  layers  over  solid  surfaces, 
the  smaller  scales  (of  the  order  of  a  grid  cell)  in  turbulent  wakes,  or  separated,  thin  vortex 
filaments  that  convect  over  arbitrarily  long  distances.  Vorticity  Confinement  can  be 
implemented  in  a  pre-existing  flow  solver,  for  both  incompressible  and  compressible 
flow,  by  adding  a  term  to  the  discretized  momentum  conservation  equations  [4],  The 
same  basic  approach,  applied  to  convecting  scalars,  will  also  be  described  below. 

For  general  unsteady  incompressible  flows,  the  governing  equations  with  the 
Vorticity  Confinement  term  are  discretizations  of  the  continuity  equation  and  the 
momentum  equations,  with  an  added  term: 


V-<?  =  0 

2.0 

d,q  =  -(q-V)q-—Vp  +  h2[y.V2q  +  ss]/At 
P 

where  q  is  the  velocity  vector,  p  is  the  pressure,  p  is  the  density,  and  p  is  a  diffusion 
coefficient  that  includes  numerical  effects  (we  assume  physical  diffusion  is  much 
smaller),  and  the  discretized  grid  cell  size  is  h  and  time  step,  At.  For  the  last  term,  £.? ,  e 
is  a  numerical  coefficient  that,  together  with  p ,  controls  the  size  and  time  scales  of  the 
convecting  vortical  regions  or  vortical  boundary  layers.  For  this  reason,  we  refer  to  the 
two  terms  in  the  brackets  as  “confinement  terms”. 

The  basic  idea  is  that  we  want  the  computed  thin  features  to  maintain  shape  and 
total  amplitude  as  their  centroid  is  convected  through  the  flow  field.  The  requirement  that 
they  relax  to  their  shape  in  a  small  number  of  time  steps  and  have  a  support  of  a  small 
number  of  grid  cells  determines  the  two  parameters,  £  and  p  .  Also,  the  “outer”  flow 
field  in  which  the  feature  is  convecting  is  slowly  varying  in  time  and  space  compared  to 
these  scales  (this  is  required  if  the  grid  cell  size  and  time  step  are  to  resolve  the  pde’s 
governing  this  outer  flow).  We  then  have  a  two-scale  problem  with  the  thin  vortex  or 
scalar  obeying  a  “fast”  dynamics. 

pV2# +es  « 0, 

where  s  is  defined  below  in  Sec.  2.1. 

Thin  vortices  are  then  convected  through  the  flow  field  by  the  “slow”  variables,  q 
(with  the  self-induced  velocity,  which  does  not  affect  the  motion  of  the  centroid  of  the 
vortex,  subtracted).  Exactly  the  same  type  of  discussion  applies  to  the  convection  of  a 
passive  scalars,  as  described  below. 

In  general,  for  convecting  vortex  filaments,  computed  flow  fields  external  to  the 
vortical  regions  are  not  sensitive  to  the  parameters  £  and  p  over  a  wide  range  of  values. 
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For  example,  the  flow  outside  an  axisymmetric  2-D  vortex  core  is  independent  of  the 
vortical  distribution,  and  hence  does  not  depend  on  e  and  p.  as  long  as  the  core  is  thin. 
Hence,  the  issues  involved  in  setting  them  are  similar  to  those  involved  in  setting 
numerical  parameters  in  other  standard  computational  fluid  dynamics  schemes,  such  as 
artificial  dissipation  in  many  conventional  compressible  solvers.  There,  the  flow  outside 
1-D  captured  shock  regions  does  not  depend  on  the  exact  internal  structure,  as  long  as  it 
is  thin. 


An  important  feature  of  the  Vorticity  Confinement  method  is  that  the 
Confinement  terms  are  non-zero  only  in  the  vortical  regions,  since  both  the  diffusion  term 
and  the  anti-diffusion  term  vanish  outside  those  regions  (care  has  to  be  taken  in  the 
numerical  implementation  to  preserve  this  feature). 

Another  important  feature  concerns  the  total  change  induced  by  the  confinement 
correction  in  mass,  vorticity  and  momentum,  integrated  over  a  cross  section  of  a 
convecting  vortex.  It  can  be  shown  [7-9]  that  mass  is  conserved  because  of  the  pressure 
projection  step  in  the  solver  and  vorticity  is  explicitly  conserved  because  of  the  vanishing 
of  the  confinement  term  outside  the  vortical  regions.  For  the  original  confinement 
formulation  used  previously,  momentum  was  almost  exactly  conserved.  A  previous 
extension  of  the  original  method  [5,10],  explicitly  conserved  the  momentum,  but  resulted 
in  a  more  spread  vortex  and  implementation  was  somewhat  complex.  The  new 
formulation,  described  below,  explicitly  conserves  momentum,  is  much  simpler,  and 
results  in  more  compact  structures. 


2.1  New  Vorticity  Confinement  Formulation 

First,  the  new  formulation  for  scalar  confinement  will  be  given.  Then,  a  velocity- 
based  (primitive  variable)  Vorticity  Confinement  term  given  that  exactly  reduces  to  the 
new  scalar  confinement  (in  terms  of  vorticity)  when  the  curl  is  taken.  The  scalar 
formulation  presented  here  is  related  to  that  presented  in  [11]  in  1-D.  In  this  section  no 
convection  is  used,  only  the  two  confinement  terms,  so  that  the  behavior  can  be  seen 
more  clearly.  Excellent  results  are  found  with  convection  and  will  be  shown  in  Sec.  3  for 
vorticity  as  well  as  convecting  scalars. 

We  start  with  an  iteration  for  a  non-negative  scalar,  <J> : 

(J)"*1  =4)n  +  p/i2V2<t>"-s/i2V2<I>"  2.1.1 


where 


<j>" 


Zcf 


♦  "H+r+s 


2.1.2 
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where  the  sum  is  over  a  set  of  grid  nodes  near  and  including  the  node  where  $  is 
computed,  the  absolute  value  is  taken  and  8  ,  a  small  positive  constant  (~  10~8)  is  added 
to  prevent  problems  due  to  finite  precision.  The  coefficients,  C,,  can  depend  on  /,  but 
good  results  are  obtained  by  simply  setting  them  all  to  1  (except  for  a  convecting  scalar, 
as  explained  in  Sec.  3.2.  Eq.  2.1.2  is  related  to  the  harmonic  mean  [13]. 

For  example,  in  2-D,  except  for  convecting  scalars,  the  form  used  in  this  study  is 


where  the  number  of  terms  in  the  sum  is  N=9.  Different  values  were  used  for  convecting 
scalars,  as  explained  in  Sec.  3.2. 

Here,  we  assume  <|>  ”  >  0 .  Negative  values  can  also  be  accommodated  with  a 
small  extension.  Both  p  and  6  are  positive. 

An  important  feature  is  that  all  terms  are  homogeneous  of  degree  1  in  Eq.  2.1.1. 
This  is  important  because  the  confinement  should  not  depend  on  the  scale  of  the  quantity 
being  confined.  Another  important  feature  is  the  nonlinearity.  It  is  easy  to  show  that  a 
linear  combination  of  terms,  for  example  of  second  and  fourth  order,  cannot  lead  to  a 
stable  confinement  for  any  finite  range  of  coefficients. 

For  smooth  <|>  fields  (long  wavelengths),  the  last  term  represents  a  diffusion  of 

order  n2  /At ,  i.e.,  it  is  first  order  accurate  if  At  is  proportional  to  h.  If  p  <  e  ,  the  total 
diffusion  (in  the  long  wavelength  limit)  is  negative.  However,  the  iteration  of  Eq.  2.1.1  is 
still  stable  and  converges  for  values  of  e  up  to  about  5  p .  Also,  the  discrete,  converged 
solution  as  «  — >  oo  can  be  given  exactly  in  terms  of  sech  functions  (to  be  described  in  the 
next  section). 

The  next  step  in  our  development  involves  letting  <|>  be  the  magnitude  of  vorticity 
and  deriving  an  equation  for  the  corresponding  primitive  variable  velocity  term  that  leads 
to  Eq.  2.1.1  when  the  curl  is  taken  (exactly  in  2-D,  or  for  a  straight  vortex  in  3-D).  This 
will  result  in  a  new  formulation  of  Vorticity  Confinement.  In  3-D,  if  the  curvature  of  a 
vortex  filament  is  large  compared  to  its  thickness,  the  results  will  still  be  approximately 
valid  since  then  the  flow  is  close  to  that  of  a  2-D  vortex  in  a  plane  normal  to  the  filament. 

We  define 

qn+'  =q"  +  \ih2V2qn  +eh2sn 

j  =  V  x  u>  2.1.3 

Since  V  •  q  =  0  ,  we  can  write 

-n+!  _  -»  _/j2Vx(p<o"  -sw") 
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where 


and 


<5"  =V*qn 


2.1.4 


where  co"  =  |co”|  +  5 

and  the  sum  is  the  same  as  in  Eq.  2. 1 .2. 

We  are  currently  implementing  the  new  method  in  our  general  fluid  dynamic 
codes.  In  addition,  the  use  of  the  “Scalar  Confinement”  version,  Eq.  2.1.1,  is  being  used 
for  flows  where  thin  streams  of  passive  scalars,  such  as  contaminants,  must  be  convected 
over  long  distances. 


2.2  Analysis  of  Zero  Convection  Form 

When  convection  and  (for  vorticity)  pressure  terms  are  included,  the  results  of 
Sec.  2.1  may  change.  For  example,  if  one  convection  step  is  followed  by  one  confinement 
step,  the  results  may  be  different,  even  in  a  frame  convecting  with  the  feature.  However, 
according  to  our  earlier  1-D  studies  [11],  we  still  expect  a  steady  state  distribution  to  be 
reached  in  the  convecting  frame.  Further,  if  desired,  we  should  be  able  to  relax  to  the 
above  zero  convection  form  each  time  step  by  performing  a  number  of  confinement 
operations  after  each  convection  operation  or  simply  take  At  small  (continuous  time, 
discrete  space  limit).  Thus,  it  is  important  to  analyze  this  form. 

Taking  the  scalar  case,  assuming  convergence  as  n->  oo ,  we  have 

V2(|0.<j)  -£<t>)  =  0, 

If  <| >  (and  hence  O )  vanishes  outside  the  feature,  we  have  p<j)0  =  sd> , 
where  the  point  (i,  j)  is  given  the  label  /  =  0 ,  or 


♦r  =« 

eN  t 


we  take 


<| sec h[ a (x,  - x0 )] sec h[ a (y;  - y0 )] ,  A,  x0,  y0  constant, 

and  where  xi  =  ih ,  y>j=  jh,h  is  the  grid  cell  size,  and  use  the  form  corresponding  to 
C,  =  1  in  Eq.  2.1.2.  We  then  have 
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<l> 


-1 


t  =AX^ 


X  _ gOhi-ax0  _j_  g-ahi+axQ 


Y  ~  gVhj-v-yo  +  ^-a/jy'+avo 


Then,  for  any  values  of  ^4,  x0  and  y0, 

sN 

(x„  +x,  +x,_lwl„+rl+r,-,)-—x,r1=o 

M' 

or 

pAA 

(e0^  +l  +  e-aA)(eaA  +l  +  e”°*) - =  0 

F 

Hence, 

1  eN  - 

ch(ah)  =  ±[(-y-l], 

2  p 

which  determines  a  .  This  implies  that  there  is  only  a  solution  for  £  >  p ,  since  N  =  9 
here.  Of  course,  stability  has  not  been  proven  here. 


2.3  Solution  with  Convection 


It  is  important  to  realize  that  the  confinement  steps  do  not  alter  the  total 
momentum,  vorticity  or  total  amount  of  a  scalar  (assuming  vorticity  or  scalar 
concentration  has  finite  support).  Also,  the  centroid  is  not  changed.  This  is  easy  to  show, 
not  only  in  the  continuum  limit  but  also  for  the  discrete  equations,  since  the  confinement 
terms  are  second  derivatives.  The  following  argument  assumes,  for  each  convection  step 
(h),  there  is  at  least  one  confinement  step  so  that  the  feature  remains  compact.  If  (|) 
represents  a  confined  passive  scalar  or  vorticity  magnitude,  then,  using  our  conservative 
convection  routine,  we  have  the  following  relationships  for  the  dynamics  of  the 
convecting  solitary  wave  (we  describe  the  2-D  case  for  simplicity): 

For  the  convecting  passive  scalar,  we  have  a  discretization  of 

d,<j)  =-V-(<7<|))  +  62V2(p<|>  -s®)/Af 
assuming  V  •  q  =  0  .  Then, 

<T’  =<t>"  -  AtVrfjrc.  •(?♦)  + *2vL04  -£<I)) 

where  discrete  operators  are  labeled. 

For  conservative  discretization,  the  total  amplitude 

<  <d  >=  2>; 

o 

is  independent  of  n.  If  we  define  the  centroid 
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<x  >■-£*, ;+;/<o> 

0 

and  the  weighted  mean  velocity 


o 


where  xtJ  is  the  (fixed)  position  vector  of  node  (i,  j),  and  4>,y  and  q(j  are  the  scalar  value 
and  the  velocity  at  that  node,  then  the  centroid  evolves  according  to: 

<X>"+I=<X>"  +A  t<q>” 

For  vortices,  the  self-induced  velocity  which  is  included  in  the  above  sum  exactly 
cancels  and,  as  in  the  passive  scalar  case,  the  qv  can  be  taken  to  be  an  externally  applied 

(irrotational)  velocity.  The  above  result  then  still  holds. 

Since  we  are,  at  this  point,  only  interested  in  the  “expectation  values”  of  scalars  or 
vorticities  for  thin  features  and  that  the  features  remain  compact,  spread  over  only  a  few 
cells,  this  Ehrenfest-type  relation  is  exactly  what  we  need.  Only  the  variables  of 
importance  are,  effectively,  solved  for.  This  shows  that  the  features,  when  isolated, 
evolve  as  particles  with  essentially  no  internal  dynamics.  However,  we  keep  the  very 
important  Eulerian  feature  that  the  number  of  features  is  not  fixed.  We  could,  for 
example,  create  additional  solitary  waves  by  inserting  a  source:  No  additional 
computational  markers  need  be  created,  as  in  Lagrangian  schemes.  For  this  study,  we 
show  that  features  can  automatically  merge  and  reduce  in  number.  This  will  be  seen  in 
the  results  of  Sec.  3.  As  for  the  earlier  Vorticity  Confinement  method,  this  property  is 
crucial  for  the  general  treatment  of  interacting  vortical  regions,  especially  in  3-D 
[5,9,10,12]. 

3.  Results 


All  results  presented  are  in  2-D.  Also,  all  convection  cases  use  a  conservative 
second-order  centered  convection  term  together  with  a  centered  diffusion  term  and  are 
first-order-explicit  in  time.  In  all  the  result  plots,  axes  are  labeled  with  grid  node  location. 
Plots  of  amplitude  are  made  using  dense  contours  with  white  grid  lines  superimposed. 


3.1  Zero  Convection 


Results  are  presented  in  Fig.  1  after  0,  8,  and  100  iterations,  for  vorticity  and 
velocity  for  a  vortex  in  2-D.  Values  used  werep  =  .2,  s  =  0  and  s  =  5p  .  In  this  figure, 
the  vorticity  contour  levels  extend  from  about  !4  of  the  maximum  initial  value  to  the 
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maximum  so  that  a  measure  of  the  size  of  the  confined  region  can  be  determined.  It  can 
be  seen  that  the  confinement  is  very  effective  and  stable. 


3.2  Passive  Scalar  Convection 

A  good  test  involves  convection  of  a  small  passive  scalar  pulse  in  solid-body 
rotating  flow  [14].  For  this  case  a  convection  step  was  followed  by  a  single  confinement 
step  (including  diffusion).  The  full  summation  of  Eq.  2.1.2  with  Ct  —  1  leads  to  a  thin 

spreading  of  the  pulses  in  the  direction  normal  to  their  motion  (This  spreading  does  not 
occur  for  vortices,  for  reasons  described  in  Sec.  3.3).  This  scalar  spreading  seemed  to  be 
due  to  the  inclusion  of  “downwind”  values,  which  is  known  to  cause  instability.  We  then 
simply  set  C,  =  0  for  downwind  points.  The  problem  was  almost  completely  cured,  and 

only  a  very  weak  spreading  remained,  which  could  easily  be  cured,  if  desired. 

Results  are  shown  in  Fig.  2a  for  a  single  pulse  convecting  through  the  1st  rotation 
(In)  and  in  Fig.  2b  through  10th  rotation.  Contours  are  plotted  from  initial  maximum 
value  to  1/3  of  that  value.  There  are  157  time  steps  between  pulse  images  (rc/8  radius)  and 
the  final  image  was  after  12,560  time  steps.  This  represented  a  travel  of  1,256  cells  or 
about  433  pulse  diameters.  The  open  circle  labels  the  starting  location.  The  time  step  was 
not  set  to  exactly  divide  the  circuit  so  we  would  not  expect  the  final  image  to  exactly 
coincide  with  the  initial  one.  However,  the  final  position  is  within  plottable  accuracy  of 
the  calculated  angle  and  radius.  The  final  radius  can  be  seen  to  be  very  close  to  the  initial 
one.  This  is  somewhat  surprising  since  the  method  was  only  first  order  in  time  and 
second  order  in  space.  We  have  not  yet  attempted  to  reduce  the  small,  residual  cross-flow 
spreading.  This  could  easily  be  accomplished  with  a  small  amount  of  additional  negative 
diffusion  in  the  cross-stream  direction.  For  the  case  computed,  we  used  p  =0.15, 

e  =  0.75 . 

Finally,  a  run  was  made  with  no  confinement  and  a  lower,  (minimal)  diffusion  set 
for  stability.  The  maximum  amplitude  decreased  to  below  the  plotting  threshold  after 
only  71/8  radians  of  travel  (position  of  the  first  image  in  Fig.  2a). 


3.3  Vorticitv  Confinement 

The  incompressible  fluid  dynamic  equations  with  Vorticity  Confinement,  Eq.  2.0, 
were  discretized  and  solved.  Each  time  step,  the  same  convection  scheme  used  for  the 
above  scalar  was  used.  A  standard  pressure  projection  method  was  used  to  enforce 
incompressibility  [15].  This  used  a  staggered  grid  for  pressure  and  velocity  and  a  box 
scheme  for  the  pressure  gradient.  The  required  Poisson  equation  for  pressure  was  solved 
using  a  direct  method  “FishPack”.  Vorticity  Confinement  terms  -  Eqs.  2.1.3  and  2.1.4, 
were  used,  with  p  =  0.05,  e=1.5p.  Results  were  very  close  to  these  for  p=0.2, 
e=1.5p,  implying  that  the  method  is  not  very  sensitive  to  the  values  of  these 
parameters.  (Further  investigaton  of  this  point  should  be  done).  The  time  step,  based  on 
grid  cell  size  and  circumferential  velocity  near  the  vortex  cores,  was  0.4  and  the  grid  had 
128x128  cells.  For  vortices,  the  full  summation  (Eq.  2.1.2)  with  all  C,  =1  was  used, 
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without  excluding  downwind  points.  This  exclusion,  useful  for  scalar  convection,  was  not 
important  here  because  the  large  circulating  velocities  near  the  vortex  cores  rapidly 
convected  any  cross-stream  spread  vorticity  around  the  cores,  keeping  them  close  to 
axisymmetric.  The  confinement  then  kept  the  cores  compact. 


3.3.1 


First,  vortices  of  opposite  sign  were  computed  with  no  external  free-stream- 
velocity,  convecting  only  under  their  own  induced  velocities.  Each  time  step  the  total 
amplitudes  and  centroids  of  both  positive  and  negative  vorticity  were  computed.  The 
velocities  induced  by  corresponding  point  vortices  at  these  locations  were  computed  on 
the  boundaries  as  velocity  boundary  conditions,  and  Dirichlet  conditions  were  used  for 
pressure. 

Results  are  presented  in  Fig.  3  for  a  sequence  of  time  steps.  The  computed  vortex 
centroids  are  convecting  at  the  same  velocity,  to  plottable  accuracy,  as  if  they  were  point 
vortices.  The  self-induced  velocity  in  each  vortex,  although  several  times  the  convecting 
velocity  of  the  centroid,  automatically  cancels  and  has  no  effect  because  the  method 
conserves  momentum,  as  explained  in  Sec.  2.3.  The  final  images  shown  result  from  1,800 
time  steps,  yet  the  vortices  are  just  as  compact  as  initially.  It  should  be  mentioned  that  if 
the  vortices  are  significantly  closer,  there  will  be  an  interaction  and  eventually  vorticity 
will  be  exchanged.  This  is  to  be  expected  since  there  are  “tails”  of  vorticity  that  (rapidly) 
decrease  with  radius  beyond  the  region  shown. 


3.3.2 


Vortices  with  the  same  sign  were  computed  as  they  rotated  around  each  other 
under  their  induced  velocities.  The  vortices  were  initially  (and  finally)  separated  by  14 
cells.  First,  the  computation  was  done  with  no  Vorticity  Confinement,  and  with  minimum 
diffusion  required  for  convective  stability.  Contours  are  displayed  in  Fig.  4  corresponding 
to  vorticity  values  between  initial  maximum  and  54  of  that  value.  It  can  be  seen  that  the 
vorticity  rapidly  diffuses.  The  computation  was  repeated  with  Vorticity  Confinement. 
Values  of  p  and  £  were  those  used  in  Sec.  3.3.1.  Vorticity  contours  between  initial 
maximum  and  54  maximum  are  displayed  for  3  different  times  in  the  1st  loop  around  each 
other  in  Fig.  4a,  and  in  the  20th  loop  in  Fig.  4b.  The  behavior  is  as  expected:  The  vortices 
can  be  seen  to  be  essentially  the  same  at  the  end  (after  7,200  time  steps)  as  initially. 

The  large  variation  in  velocity  in  the  space  between  the  vortices  can  be  seen  in 
Fig.  5,  both  after  1  and  after  20  loops.  This  would  make  it  difficult  for  a  pde-based  Taylor 
expansion  to  give  accurate  results  over  long  convection  times. 

As  stated,  the  main  reason  for  using  an  Eulerian  method  as  opposed  to  Lagrangian 
is  to  allow  vortices  to  merge  and  otherwise  change  topology,  so  that  complex,  realistic 
vortical  configurations  can  be  treated.  This  has  been  amply  demonstrated  by  the  original 
Vorticity  Confinement  method  both  in  2-D  and  3-D  [5,7,9,10,12].  To  demonstrate  that 
the  new  confinement  method  does  not  inhibit  this  merging  we  repeat  the  above 
computation,  but  with  the  vortices  initially  5  cells  apart.  This  is  a  repeat  of  a  calculation 
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using  the  original  Vorticity  Confinement  method  [7].  It  is  known  from  vorticity  contour 
or  “waterbag”  methods  that  corotating  vortices  become  unstable  and  merge  when  they  are 
initially  closer  than  a  few  diameters.  It  can  be  seen  in  the  vorticity  contour  plots  of  Fig.  6 
that  this  merging  automatically  occurs,  as  expected,  with  the  new  Vorticity  Confinement 
method. 


Conclusion 


A  new  Eulerian  Vorticity  Confinement  has  been  introduced.  Its  advantages  over 
the  original  Vorticity  Confinement  method  are  that  it  is  simpler  and  easier  to  analyse  and 
explicitly  conserves  total  momentum.  The  method  accurately  convects  thin  features, 
either  vortices  or  scalars,  over  very  long  distances,  even  though  they  are  confined  to  only 
2-3  grid  cells.  As  in  shock  and  other  discontinuity  capturing  schemes,  the  details  of  the 
internal  structure  of  the  feature  are  not  solved  for,  but  implicitly  modeled  (the  features  are 
treated  as  weak  solutions).  Integral  quantities,  such  as  total  amplitude  and  centroid, 
however,  are  accurately  solved  for.  As  such,  the  method  essentially  represents  the 
features  as  nonlinear  discrete  solitary  waves  that  live  on  the  computational  lattice.  It  is 
argued  that  the  new  method  is  a  rational  generalization  of  1-D  discontinuity  confinement 
schemes  to  multiple  dimensions. 

Examples  are  presented  for  thin,  convecting  passive  scalars  and  vortices  in  2-D 
incompressible  flow. 
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(a)  1st  loop  (b)  10th  loop 


Figure  2  Scalar  convection 


time  step=0.4/  time  steps  between  two  plots=300 


Figure  3  Vorticity  contour  plots  of  self-induced  flow  by  two  point  vortices(15  cells 
apart)  with  opposite  sign 
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(b)  (c) 

Figure  4  Vorticity  contour  plots  of  self-induced  flow  by  two  vortices(14  cells  apart) 
with  same  sign,  (a)  without  Vorticity  Confinement;  (b)with  Vorticity  Confinement,  1 
loop;  (c)with  Vorticity  Confinement,  after  20th  loops. 


Figure  5  Vector  plots  of  self-induced  flow  by  two  vortices  (14  cells  apart)  with  same 
sign,  (a)  with  Vorticity  Confinement,  after  1  loop;  (b)  with  Vorticity  Confinement,  after 
20  loops. 
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Appendix 


The  new  Vorticity  Confinement  method  can  be  shown  to  be  related  to,  and  a 
generalization  of,  some  recent  1-D  discontinuity  preserving  schemes:  In  this  Appendix 
we  show  that  (a)  For  incompressible  flow,  the  proper  generalization  to  multidimensions 
of  the  first  order  derivatives  used  in  1-D  discontinuity  preserving  schemes  is  the  vorticity, 
which  is  the  only  non-zero  rotationally  invariant  first  order  velocity  derivative  in  2-D; 
and  (b),  The  proper  generalization  to  multidimensions  of  the  discontinuous  logical  “min” 
functions  used  by  1-D  schemes  is  the  smooth  algebraic  form,  Eq.  2.1.4,  used  in  the  new 
Vorticity  Confinement  method. 

1.  An  algebraic  nonlinear  negative  diffusion  term  similar  to  that  used  here 
was  added  to  a  first-order  upwind  method  (which  has  an  implicit  positive  diffusion)  for  1- 
D  passive  scalar  convection  [11].  The  result  for  a  constant  speed  pulse  was  a  convecting, 
stable  solitary  wave,  spread  over  a  few  grid  cells.  A  sequence  of  amplitudes  was  plotted 
for  a  large  number  of  time  steps  in  a  frame  moving  with  the  pulse,  and  it  was  shown  that 
there  was  a  smooth  underlying  structure:  The  sequence  of  amplitudes  at  the  grid  nodes  at 
different  time  steps  was  equivalent  to  that  generated  by  moving  the  smooth  structure 
through  the  grid  at  constant  speed  and  sampling  the  values  on  the  grid  nodes  at  each  time 
step.  The  current  method  appears  to  give  similar  results. 

2.  A  more  compact  method  was  derived  also  in  [11]  where  the  pulse  had 
support  of  exactly  two  grid  nodes  (two  non-zero  amplitude  values)  at  each  time  step.  This 
was  actually  derived  by  starting  from  a  Lagrangian  particle  representation  with  which  it 
is  in  1-1  correspondence,  although  it  was,  of  course,  Eulerian. 

3.  In  the  same  paper,  it  was  shown  that  when  the  scheme  was  formulated  in 

terms  of  a  new  variable, 

/ 

W/=ZW7 

0 

where  u,  was  the  original  variables.  Step  functions  could  be  propagated  indefinitely. 

There  was  a  1-1  correspondence  between  the  two  cases. 

4.  The  quantities  used  in  these  pulse  schemes  were  first  order  derivatives  of 
the  variables  used  in  the  step  function  schemes.  An  analog  of  Vorticity  Confinement  can 
be  formulated  in  these  1-D  studies  by  considering  the  summed  variable  (w)  to  be  a 
velocity  (v)  normal  to  the  1-D  axis  (x)  and  the  vorticity  the  derivative  of  this  along  the 
axis: 

co,  =  5,v(x) 

We  then  immediately  have  a  way  to  generalize  to  multi-dimensions.  Instead  of 
concatenating  1-D  operators,  we  realize  that  for  example,  in  2-D  incompressible  flow 
there  are  only  2  rotationally  invariant  first  order  derivatives  of  velocity: 

D  =  V  q  , 

®=Vx$|* 

where  k  denotes  the  direction  normal  to  the  2-D  plane. 
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For  incompressible  flow,  D  is,  of  course,  zero  and  we  are  led  to  formulate  our 
multidimensional  method  only  in  terms  of  oo.  This  is  a  rational  generalization  of  the  1-D 
5xv  term. 

5.  The  logical  “min”  function  used  in  the  compact  1-D  scheme  does  not  seem 
to  be  suitable  for  multidimensions.  While  not  important  in  1-D,  the  field  within  a  2-D  or 
3-D  feature  should  be  relatively  smooth  and  free  of  large  variations  caused  by 
discontinuous  logical  functions  as  the  feature  is  rotated.  A  final  generalization,  which 
leads  us  to  our  new  Vorticity  Confinement  method,  involves  the  function. 

/ 

If  p  is  large,  d>  -»  min^} .  We  use  p= 1  in  our  method  (Eq.  2.1.2),  which  gives  smooth, 

algebraic  results  and  approximates  the  “min”  function. 

6.  The  discontinuity  capturing  version  of  the  above  1-D  compact  pulse 
scheme  (with  “min”  functions)  of  Ref.  [1 1]  is  apparently  similar  to  the  recent  “superbee” 
scheme  [13]. 
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Vorticity  Confinement  -  Recent 
Results:  Turbulent  Wake 
Simulations  and  A  New, 
Conservative  Formulation 

J.  Steinhoff  M.  Fan  2,  and  L.  Wang  3 


Abstract 


Vorticity  Confinement  has  been  shown,  over  the  last  several  years,  to  be  a 
way  to  quickly  and  cheaply  approximate  incompressible  flows  over  complex 
configurations  at  high  Reynolds  number.  For  these  flows,  there  are  many 
important  salient  features  that  are  reproduced  by  the  method  in  a  simple  way, 
on  relatively  coarse  uniform  Eulerian  computational  grids  with  fast  low-order 
computational  methods.  No  complex,  high  order  schemes  and  no  extensive 
refinement  or  body  conforming  grid  are  required.  The  basic  features  of  the  flow 
that  are  very  effectively  simulated  all  involve  vortical  regions:  thin  boundary 
layers  on  solid  surfaces  that  are  attached  or  separate,  thin  vortex  sheets  and 
thin  filaments  that  can  convect  over  long  distances  with  no  significant  physical 
diffusion,  and,  finally,  turbulent  wakes,  including  small  scale  effects.  These 
features  are  very  difficult  to  treat  with  conventional  schemes.  The  main  point 
is  that  almost  all  of  the  vortical  regions  in  these  flows  are  either  thin,  so 
that  their  internal  structure  is  not  important  (analogous  to  shocks  or  contact 
discontinuities  in  compressible  flow),  or  contain  small  scale  structures,  such 
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as  turbulent  wakes.  All  of  these  vortical  structures  are,  of  course,  embedded 
in  an  incompressible  irrotational  flow  field,  which  is  also  efficiently  computed. 

First,  the  salient  features  of  Vorticity  Confinement  solutions  will  be 
discussed,  to  give  the  reader  an  understanding  of  the  characteristics  of  the 
method.  Then,  a  short  presentation  of  the  original  formulation  of  the  method 
will  be  given  and  a  number  of  recent  papers  giving  more  details  referenced. 
This  will  be  followed  by  some  recent  results  where  the  method  is  used  to  treat 
the  small  scales  in  turbulent  wakes  in  LES-type  simulations.  Many  other  recent 
results  concern  free  vortices  and  attached  flow  over  complex  configurations, 
as  well  as  validation  studies.  Papers  describing  these  will  be  referred  to. 
Finally,  a  new,  simpler  formulation  will  be  described  that  is  effective  for  both 
Vorticity  Confinement  as  well  as  convection  of  thin  streams  of  passive  scalars. 
Preliminary  results  of  this  new  formulation  for  convecting  vortices  and  scalars 
in  2-D  will  then  be  presented. 


1.1  Introduction 

Vorticity  Confinement  has  been  shown,  over  the  last  several  years,  to  be  a 
way  to  quickly  and  cheaply  approximate  incompressible  flows  over  complex 
configurations  at  high  Reynolds  number^1  ~10).  For  these  flows,  there  are  many 
important  salient  features  that  are  reproduced  by  the  method  in  a  simple  way, 
on  relatively  coarse  uniform  Eulerian  computational  grids  with  fast  low-order 
computational  methods.  No  complex,  high  order  schemes  and  no  extensive 
refinement  or  body  conforming  grid  are  required.  The  basic  features  of  the  flow 
that  are  very  effectively  simulated  all  involve  vortical  regions:  thin  boundary 
layers  on  solid  surfaces  that  are  attached  or  separate,  thin  vortex  sheets  and 
thin  filaments  that  can  convect  over  long  distances  with  no  significant  physical 
diffusion,  and,  finally,  turbulent  wakes,  including  small  scale  effects.  These 
features  are  very  difficult  to  treat  with  conventional  schemes.  The  main  point 
is  that  almost  all  of  the  vortical  regions  in  these  flows  are  either  thin,  so 
that  their  internal  structure  is  not  important  (analogous  to  shocks  or  contact 
discontinuities ( 1 1 )  in  compressible  flow),  or  contain  small  scale  structures,  such 
as  turbulent  wakes.  All  of  these  vortical  structures  are,  of  course,  embedded 
in  an  incompressible  irrotational  flow  field,  which  is  also  efficiently  computed. 

The  main  idea  behind  Vorticity  Confinement  is  that,  at  high  Reynolds 
number,  the  above  small  vortical  scales  are  modeled  on  the  grid  in  an 
inherently  discrete,  simple  way  that  is  much  more  efficient  than  if  model 
partial  differential  equations  (pde’s)  were  first  formulated  and  then  resolved 
with  finite  difference  approximations.  This  modeling  involves  nonlinear  terms 
in  the  discrete  equations  that  result  in  a  stable  negative  total  diffusion  for  the 
vorticity  in  a  band  of  length  scales.  This  negative  diffusion  acts  to  confine  the 
small  scale  vorticity  to  a  thickness  of  about  2  grid  cells:  Below  this  scale, 
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the  structure  is  dominated  by  a  positive  diffusion.  On  scales  significantly 
larger,  the  negative  diffusion  causes  the  structure  to  contract.  The  equilibrium 
solution  results  from  a  balance  of  these  two  effects.  An  important  feature  is 
that  the  irrotational  flow  regions,  as  well  as  the  very  large  vortical  scales 
that  are  present  in  the  turbulent  wakes  are  not  significantly  affected  by  these 
Vorticity  Confinement  terms  and  are  accurately  resolved  and  convected  by  the 
method,  which  reduces  to  a  conventional  CFD  method  at  these  large  scales. 

First,  the  salient  features  of  Vorticity  Confinement  solutions  will  be 
discussed,  to  give  the  reader  an  understanding  of  the  characteristics  of  the 
method.  Then,  a  short  presentation  of  the  original  formulation  of  the  method 
will  be  given  and  a  number  of  recent  papers  giving  more  details  referenced. 
This  will  be  followed  by  some  recent  results  where  the  method  is  used  to 
treat  the  small  scales  in  turbulent  wakes  in  LES-type  simulations^.  Many 
other  recent  results  concern  free  vortices  and  attached  flow  over  complex 
configurations,  as  well  as  validation  studies.  Papers  describing  these  will  be 
referred  to.  Finally,  a  new,  simpler  formulation  will  be  described  that  is 
effective  for  both  Vorticity  Confinement  as  well  as  convection  of  thin  streams 
of  passive  scalars.  Preliminary  results  of  this  new  formulation  for  convecting 
vortices  and  scalars  in  2-D  will  then  be  presented^ . 


1.2  Vorticity  Confinement:  Salient  Features 

There  are  several  different  ways  to  characterize  Vorticity  Confinement  -  they 
are  all  important  in  understanding  how  it  works: 

1.  It  is  a  way  of  treating  flow  with  small  vortical  scales  as  “weak  solutions” 
of  the  Euler  equations,  where  the  small  scales  are  treated  as  regularized 
singular  regions.  In  this  way  it  is  analogous  to  shock  capturing-quantities  are 
conserved,  not  at  each  point  but  integrated  through  the  vortical  feature,  which 
is  spread  over  a  few  grid  cells.  The  confinement  terms  serve  as  a  simple  implicit 
model  of  the  internal  structure  and  have  little  effect  on  the  flow  external  to 
the  feature,  which  is  accurately  treated.  Well  known  examples  of  this  type 
of  treatment  include,  besides  shock  capturing,  vortex  filament  and  vortex 
lattice  schemes  for  convecting  vortex  sheets.  Also,  the  use  of  panel  methods 
for  external  flows,  where  thin  vortical  boundary  layers  are  simply  treated 
involves  similar  ideas(12-16).  An  important  example  of  this  independence  of 
the  flow  outside  the  features  involves  vortex  filaments  with  near-axisymmetric 
cores  with  radius  of  curvature  much  larger  than  core  radius,  where  the  flow  is 
approximately  2-D  in  planes  normal  to  the  filaments. 

2.  The  current  Vorticity  Confinement  method  can  be  treated  as  a  “zeroth 
order”  approximation  for  the  thin  vortical  features.  As  in  well-known 
interacting  boundary  layer  schemes,  the  basic  solution  can  be  treated  as 
a  first  term  in  an  asymptotic  expansion.  Higher  order  terms  would  involve 
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modeling  some  of  the  effects  of  the  structure,  such  as  displacement  thickness 
and  skin  friction.  In  addition  to  boundary  layers,  such  secondary  effects  could 
include  long-term  spreading  of  aircraft  trailing  vortices<17)  and  modeling  of 
any  effects  of  axial  flow  on  vortex  breakdown.  For  most  of  the  problems  that 
we  have  treated  until  now  with  Vorticity  Confinement,  these  effects  have  not 
been  very  important  and  the  “zeroth  order”  approach  was  sufficient.  Cases 
where  such  higher  order  effects  have  already  been  treated  include  long-term 
trailing  vortex  simulation^  17\  and  a  study  involving  flow  separating  from  a 
curved  surface^  where  the  boundary  layer  was  turbulent  and  affected  the 
separation  location. 

It  should  be  mentioned  that,  of  course,  a  very  large  amount  of  work  has 
been  done,  over  the  last  few  decades,  on  turbulent  boundary  layer  modeling 
using  conventional  pde-based  models,  (mostly  calibrating  Reynolds  averaged 
“RANS”  models).  The  same  is  true  of  the  small  scales  in  turbulent  wakes, 
using  large  eddy  simulation  (LES).  By  comparison,  very  little  turbulent 
modeling  has  been  done  based  on  Vorticity  Confinement,  i.e.,  by  adding 
“higher  order”  additional  discrete  terms  with  adjustable  parameters  to  the 
basic  method.  For  these  reasons  we  cannot  claim  that  Vorticity  Confinement 
turbulence  models  currently  exist  for  these  “first  order”  effects.  We  are  only 
stating  that  Vorticity  Confinement  appears  to  be  a  promising  approach  and 
a  framework  for  modeling,  because  of  its  simplicity  and  efficiency.  Additional 
work  should  be  done  to  develop  confinement-based  turbulence  models. 

3.  For  convecting  vortex  filaments,  the  internal  dynamics  is,  effectively, 
treated  as  “fast”  dynamics  which  approximately  satisfies  the  confinement 
terms  separately  from  the  other,  conventional  CFD  terms  used  in  the  method 
(these  will  be  described  below).  This  fast  dynamics  is  then  slaved  to  the  basic 
fluid  dynamics  as  the  filaments  are  “slowly”  convected  through  the  flow. 

4.  The  Vorticity  Confinement  method  can  be  thought  of  as  incorporating 
both  conventional  CFD  ideas  for  solving  for  the  “smooth”,  large  scale  features 
of  the  flow,  and  intrinsically  discrete  ideas  for  solving  for  the  small  scale 
features:  The  confinement  terms  result  in  intrinsically  discrete  thin,  small 
scale  solitary-wave  structures  that  “live”  on  the  computational  lattice.  In  this 
way,  it  is  similar  to  lattice  gas  approaches.  For  the  large  scales,  as  mentioned, 
the  confinement  terms  have  little  effect  and  the  method  reverts  to  conventional 
CFD,  which  is  well-known  to  be  effective  for  these  scales. 

For  small  scales,  the  method  corrects  the  well  known  difficulty  of  treating 
thin  features  (with,  of  course,  large  gradients)  inherent  in  conventional  CFD 
schemes.  For  example,  very  fine  body  conforming  grids  are  required  with 
conventional  “Navier-Stokes”  turbulence  modeling  methods  for  boundary 
layers,  and  even  then,  the  result  is  a  pde-based  model  solution.  Thin  boundary 
layers  are  very  simply  treated  with  confinement  as  intrinsically  discrete 
structures. 

On  the  other  hand,  while  intrinsically  discrete  models  such  as  “lattice  gas” 
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can,  perhaps,  be  efficiently  formulated  for  thin  regions  such  as  turbulent 
boundary  layers,  their  treatment  of  large  scales  is  very  inefficient.  For  example, 
lattice  gas  schemes  only  converge  statistically  to  smooth  laminar  Navier-stokes 
solutions;  Samples  of  the  solution  at  large  numbers  of  grid  points  (iV5),  over 
large  numbers  of  time  steps  (Nt)  must  be  used  in  any  small  region  in  time  and 
space  in  order  to  estimate  the  solution  there.  For  large  scales,  therefore,  these 
methods  converge  like  s/NsNt  and  can  be  thought  of  as  “half-order”  accurate, 
as  opposed  to  conventional  CFD,  for  which  very  high  order  methods  have  been 
developed^18*. 


1.3  Vorticity  Confinement:  Basic  Formulation 

A  specific  solver  will  be  described  which  employs  Vorticity  Confinement. 
However,  Vorticity  Confinement  can  be  implemented  in  a  pre-existing  flow 
solver,  for  both  incompressible  and  compressible  flow,  by  adding  a  term  to 
the  discretized  momentum  conservation  equations^. 

For  general  unsteady  incompressible  flows,  the  governing  equations  with  the 
Vorticity  Confinement  term  are  discretizations  of  the  continuity  equation  and 
the  momentum  equations,  with  an  added  term: 


V-q  =  0 


3.0 


dtq  —  —(q-V)q  +  V(p/p)  +  [pV2q  —  €s\h2/At 


where  q  is  the  velocity  vector,  p  is  the  pressure,  p  is  the  density,  and  p  is 
a  diffusion  coefficient  that  includes  numerical  effects  (we  assume  physical 
diffusion  is  much  smaller).  For  the  last  term,  es,  e  is  a  numerical  coefficient 
that,  together  with  p,  controls  the  size  and  time  scales  of  the  convecting 
vortical  regions  or  vortical  boundary  layers.  For  this  reason,  we  refer  to  the 
two  terms  in  the  brackets  as  “confinement  terms”.  Also,  h  is  the  grid  cell  size 
and  At  the  time  step. 

The  first  confinement  term,  or  diffusion  term  is  usually  specified  explicitly, 
but  can  be  implicitly  present  in  the  solver,  such  as  in  a  lower  order  upwind 
scheme.  There  are  also  many  possible  forms  for  the  second  confinement  term. 

First,  the  original  one  used  in  this  and  earlier  studies  will  be  described.  Then, 
a  new  more  elegant,  simpler  form  will  be  described,  together  with  a  simple 
demonstration. 
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1.4  Original  Vorticity  Confinement 

1.4.1  Formulation 

We  define  the  second  confinement  term  as 

_  - 
s  =  —n  x  uj 
h 

where  _ 

.  _  Vq 

n  ~  |V»?| 

and  the  vorticity  vector  is  given  by 

u  =  V  x  q 


and 

V  =  |w| 

In  general,  for  boundary  layers  and  convecting  vortex  filaments,  since 
computed  flow  fields  external  to  the  vortical  regions  are  not  very  sensitive 
to  the  details  of  the  internal  structure  as  long  as  they  are  smooth,  they 
are  not  sensitive  to  the  parameters  e  and  /x  over  a  wide  range  of  values. 
Hence,  the  issues  involved  in  setting  them  are  similar  to  those  involved  in 
setting  numerical  parameters  in  other  analogous  computational  fluid  dynamics 
schemes,  such  as  artificial  dissipation  in  conventional  shock  capturing  schemes. 
For  wake  flows,  e  and  /x  can  be  used  to  approximately  simulate  finite  Reynolds 
number  effects,  since  they  control  the  intensity  of  the  small  vortical  scales. 
This  is  an  area  of  study  that  is  just  beginning  -  some  results  of  which  will  be 
reported  here. 

An  important  feature  of  the  Vorticity  Confinement  method  is  that  the 
confinement  terms  are  non-zero  only  in  the  vortical  regions,  since  both  the 
first  (diffusion)  term  and  the  second  (anti-diffusion)  term  vanish  outside  those 
regions  (care  has  to  be  taken  in  the  numerical  implementation  to  preserve  this 
feature). 

Another  important  feature  concerns  the  total  change  induced  by  the 
confinement  correction  in  mass,  vorticity  and  momentum,  integrated  over  a 
cross  section  of  a  convecting  vortex.  It  can  be  shown^5,9^  that  mass  is  conserved 
because  of  the  pressure  projection  step  in  the  (incompressible)  solver  and 
vorticity  is  explicitly  conserved  because  of  the  vanishing  of  the  correction 
outside  the  vortical  regions.  Momentum  is  almost  exactly  conserved.  An 
extension  of  the  original  method^ 10)  explicitly  conserves  the  momentum  at 
the  expense  of  some  additional  complications.  The  new,  simpler  formulation, 
described  below  in  Sec.5.2  explicitly  conserves  momentum  with  no  required 
extensions. 
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1.4.2  Results  of  Wake  Study 

Two  test  cases  are  presented  in  this  section  involving  the  wake  structure  of  a 
3-D  circular  and  square  cylinder. 


1.4. 2.1  3-D  Circular  Cylinder 

Flow  over  a  3-D  circular  cylinder  was  calculated  to  assess  the  ability  of 
Vorticity  Confinement  to  accurately  model  the  wake  flow  behind  a  blunt 
body.  A  long  cylinder  was  "immersed”  in  a  uniform  141  x  101  x  61  Cartesian 
grid  in  the  streamwise,  normal  and  spanwise  directions,  respectively.  Periodic 
conditions  were  imposed  at  the  lateral  boundaries.  The  diameter  of  the 
cylinder  was  15  grid  cells.  Comparisons  between  experiment  and  computed 
results  are  given  in  Figs.  1  and  2.  In  the  figures,  the  origin  of  the  coordinate 
system  used  is  located  in  the  center  of  the  cylinder,  and  all  distances  are  non- 
dimensionalized  by  the  diameter.  The  experimental  results  of  Lourenco  and 
Shih<19>  at  a  Reynolds  number  of  about  3,900,  were  compared  to. 

The  diffusion  coefficient  \i  was  held  constant  for  this  study.  The  confinement 
coefficient,  e,  was  adjusted  to  impose  different  levels  of  confinement.  This 
resulted  in  different  levels  of  the  intensity  of  the  small  vortical  scales  in  the 
wake,  and  approximately  simulated  different  Reynolds  numbers. 

Figure  1  depicts  the  mean  streamwise  velocities  resulting  from  three  values 
of  the  confinement  coefficient.  Figures  a)  and  b)  depict  the  result  of  two  levels 
of  confinement,  where  the  flow  is  measured  along  lines  traversing  the  wake 
(normal  to  the  cylinder  axis  and  the  mean  stream)  at  three  different  locations 
in  the  wake  of  the  cylinder.  For  both  values  of  e  the  agreement  can  be  seen 
to  be  very  good,  indicating  that  the  effect  of  the  confinement  parameter  is 
small  over  a  range  of  values.  Results  from  a  case  without  confinement  (e  =  0) 
are  depicted  in  Fig.  2c.  Without  confinement,  the  flow  field  is  dominated  by 
diffusive  effects  that  are  not  counterbalanced  by  the  anti-diffusive  confinement 
term  and  approximates  a  steady,  low  Reynolds  number  flow. 

The  ability  of  Vorticity  Confinement  to  model  the  turbulent  wake  was 
assessed  by  computing  the  rms  streamwise  velocity  fluctuations  in  the  wake 
region.  Comparisons  of  these  fluctuations  with  the  experimental  data  were 
made  along  the  same  lines  in  the  wake  as  for  the  mean  velocity,  for  the  same 
three  values  of  e.  Fig.  2a  (e  =  0.25)  shows  very  good  agreement.  This  is 
the  same  value  that  shows  the  best  agreement  for  the  mean  velocity,  as  can 
be  seen  in  Fig.  1.  The  effect  of  increasing  confinement  to  0.5  is  depicted  in 
Fig.  2b.  In  general,  the  effect  of  increased  confinement  is  to  thin  the  shear 
layer  comprising  the  wake  boundary,  and  to  increase  the  fluctuation  of  the 
time-dependent  flow  from  the  mean  flow.  Fig.  3c  depicts  results  without 
confinement.  Without  confinement,  the  flow  can  be  seen  to  be  steady  and 
exhibits  none  of  the  fluctuations  that  occur  when  confinement  is  used. 
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Figure  3  depicts  isosurfaces  of  vorticity  magnitude  for  the  same  three  levels 
of  confinement.  The  use  of  confinement  results  in  chaotic  flow  patterns,  as 
would  be  expected  in  three-dimensional  turbulent  flows.  This  does  not  occur  in 
two-dimensional  simulations,  indicating  that  this  chaotic  behavior  is  not  due 
to  numerical  instability  created  by  the  confinement.  Increasing  confinement 
increases  the  chaotic  nature  of  the  flow  in  3-D  and  reduces  the  characteristic 
size  of  the  vortical  structures,  analogous  to  what  would  be  expected  with  an 
increase  in  Reynolds  number.  Clearly,  the  use  of  confinement  allows  small- 
scale  time-dependent  wake  structures  to  be  generated  on  extremely  coarse 
grids  with  the  small-scale  structure  captured  over  only  1  ~  2  grid  cells.  Also, 
these  small-scale  structures  serve  as  a  viscous  sink  for  turbulent  energy,  as  in 
physical  turbulence. 


1. 4.2.2  3-D  Square  Cylinder 

Flow  over  a  square  cylinder  was  also  calculated.  As  in  the  circular  case,  the 
cylinder  was  “immersed”  in  a  uniform  141  x  101  x  61  Cartesian  grid  and 
periodic  conditions  imposed  at  the  lateral  boundaries.  The  diameter  (length 
of  each  side)  of  the  cylinder  was  also  15  grid  cells.  The  same  coordinate 
system  was  used  as  for  the  circular  cylinder.  Results  of  the  computations 
were  compared  to  the  experimental  results  of  Lyn  et  al/20^  at  a  Reynolds 
number  of  about  21,400. 

As  in  the  circular  cylinder  case,  the  diffusion  coefficient  (i  was  set  to  0.15. 
The  confinement  coefficient,  e,  was  adjusted  to  impose  different  levels  of 
confinement  so  as  to  approximate  the  effects  of  different  Reynolds  numbers. 

Figure  4  depicts  the  comparison  with  experimental  data  of  the  time- 
averaged  streamwise  velocity  along  a  streamwise  line  extending  downstream 
from  the  middle  of  the  leeward  face  of  the  cylinder.  Results  of  two  values  of  the 
confinement  coefficient  are  plotted.  Figure  6  shows  the  time-averaged  velocity 
along  a  line  traversing  the  wake  (normal  to  the  cylinder  axis  and  the  mean 
stream),  as  in  the  circular  case,  at  x  =  1.  Symbols  represent  the  experimental 
data.  Our  numerical  results  agree  well  with  the  experimental  data  for  e  =  0.35 
and  e  =  0.25. 

Comparisons  of  the  computed  turbulence  level  with  the  experimental  results 
also  show  reasonably  good  agreement  in  Fig.  6a.  As  in  the  circular  case, 
the  same  value  of  e  that  gave  good  agreement  with  experiment  for  RMS 
velocity  also  gave  good  agreement  formean  velocity.  The  effect  of  decreasing 
confinement  is  shown  in  Fig.  6b.  Comparing  this  case  with  the  previous  3-D 
circular  cylinder  case,  it  is  easy  to  see  that  a  higher  value  of  e  is  required  for 
better  results.  This  is  apparently  because  the  Reynolds  number  has  increased 
from  about  3,900  in  the  3-D  circular  cylinder  case  to  roughly  21,400  in  the 
present  3-D  square  cylinder  case. 
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I.4. 2. S  Synopsis  of  Cylinder  Study 

In  subsequent  studies,  the  correlation  between  e  and  Reynolds  number  will  be 
studied  by  comparisons  between  computation  and  experiment  for  other  cases 
and  a  useful  calibration  of  e  determined.  Also,  it  should  be  emphasized  that 
our  agreement  with  experiment  is  closer  than  many  much  finer,  conforming 
grid  studies  with  much  more  complex  LES  pde  models^21),  each  requiring  a 
number  of  empirical  coefficients.  As  explained,  just  one  coefficient  is  adjusted 
here,  on  a  simple,  coarse  uniform  Cartesian  grid. 


1.5  New  Vorticity  Confinement 

1.5.1  Formulation 

Because  additional  corrections  must  be  added  to  make  the  original 
confinement  explicitly  conserve  momentum,  a  new,  simpler  formulation  that 
does  not  have  these  problems  has  been  developed.  This  new,  intrinsically 
discrete  formulation  is  presented  in  this  section. 

A  more  detailed  description  is  presented  in  [1].  First,  a  formulation  for  scalar 
confinement  is  given.  Then,  a  velocity-based  (primitive  variable)  Vorticity 
Confinement  correction  is  given  that  reduces  to  the  scalar  confinement  in 
terms  of  vorticity  when  the  curl  is  taken. 

1.5. 1.1  Scalar  Formulation 

The  scalar  formulation  presented  here  is  related  to  that  presented  in  [8]  in 
1-D. 

We  start  with  the  form  for  a  scalar  : 

<j>n+l  =  4>n  -  A tV  •  #n  +  6<j)n  5.1.1 

where  the  confinement  correction: 

6<f>n  =  /i2[/xVVn  ~  eV2$n] 


where 


r£  ciir+sr1] 


-1 


5.1.2 


where  the  sum  is  over  a  set  of  grid  nodes  near  and  including  the  node  where  $ 
is  computed,  and  <5  is  a  small  positive  constant  (~  10  8)  to  prevent  problems 
due  to  finite  precision.  The  coefficients,  Ci,  can  be  varied  but  good  results  are 
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obtained  by  simply  setting  them  to  1.  For  example,  in  2-D,  one  possibility  is 

a=-l  /3=-l _  r  1  Q 


where  N  neighboring  terms  are  taken. 

Here,  we  assume  (j>n  >  0.  Negative  values  can  also  be  accommodated  with 
a  small  extension.  Both  fi  and  e  are  positive. 

An  important  feature  is  that  all  terms  are  homogeneous  of  degree  1  in  Eq. 
5.1.1.  This  is  important  because  the  confinement  not  depend  on  the  scale  of 
the  quantity  being  confined.  Another  important  feature  is  the  nonlinearity. 
It  is  easy  to  show  that  a  linear  combination  of  terms  of  different  order  in 
the  derivatives  cannot  lead  to  a  stable  confinement  for  any  finite  range  of 
coefficients. 

For  smooth  <j>  fields  (large  scales),  the  last  term  represents  a  negative 
diffusion.  If  jj  <  e,  where  N  is  the  number  of  terms  in  Eq.  5.1.2,  the  total 
diffusion  is  negative.  However,  the  iteration  of  Eq.  5.1.1  is  still  stable  and 
converges  for  values  of  e  up  to  about  5/x,  resulting  in  an  effective  negative 
diffusion  coefficient  in  the  long  wavelength  limit.  Also,  for  q  =  0,  the  discrete, 
converged  solution  can  be  given  exactly  in  terms  of  (translationally  invariant) 
sech  functions. 

For  convection  using  a  conservative  discretization,  the  total  amplitude 

ij 

is  independent  of  n.  If  we  define  the  centroid  of  an  isolated  region, 

<x>n=j2^j/  <^> 

ij 


and  the  weighted  mean  velocity 

<  q  >"=  <  *  > 

ij 

where  Xij  is  the  (fixed)  position  vector  of  node  (i,  j),  and  faj  and  are  the 
scalar  value  and  the  velocity  at  that  node,  then  the  centroid  evolves  according 
to: 

<  x  >n+1=<  X  >n  +A t  <Q>n 

For  vortices,  the  self-induced  velocity  which  is  included  in  the  above  sum 
exactly  cancels  and,  as  in  the  passive  scalar  case,  the  qfij  can  be  taken  to  be 
an  externally  applied  (irrotational)  velocity.  The  above  result  then  still  holds. 
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Since  we  are,  at  this  point,  only  interested  in  the  “expectation  values”  of 
scalars  or  vorticities  for  thin  features  and  that  the  features  remain  compact, 
spread  over  only  a  few  cells,  this  Ehrenfest-type  relation  is  exactly  what  we 
need.  Only  the  variables  of  importance  are,  effectively,  solved  for.  This  shows 
that  the  features,  when  isolated,  evolve  as  particles  with  essentially  no  internal 
dynamics.  However,  we  keep  the  very  important  Eulerian  feature  that  the 
number  of  features  is  not  fixed.  We  could,  for  example,  create  additional 
solitary  waves  by  inserting  a  source:  No  additional  computational  markers 
need  be  created,  as  in  Lagrangian  schemes.  For  this  study,  we  show  that 
features  can  automatically  merge  and  reduce  in  number.  This  will  be  seen  in 
the  results  of  Sec.  5.2.  As  for  the  earlier  Vorticity  Confinement  method,  this 
property  is  crucial  for  the  general  treatment  of  interacting  vortical  regions, 
especially  in  3-D(4,7,9,10). 

1.5. 1.2  Vorticity  Formulation 

The  next  step  involves  letting  <j)  be  the  magnitude  of  vorticity  and  deriving  an 
equation  for  the  corresponding  velocity  correction  that  leads  to  Eq.  5.1.1  when 
the  curl  is  taken  (in  2-D).  This  correspondence  is  still  close  for  a  3-D  vortex 
filament  if  the  radius  of  curvature  is  large  so  that  the  flow  is  approximately  2- 
D  in  sections  normal  to  the  filament  axis.  This  will  result  in  a  new  formulation 
of  Vorticity  Confinement. 

We  simply  define 

^  -F  AtV  *  <f  qn  +  Sq"  5.1.4 

where  the  confinement  correction: 

5 q"  =  h 2(fiV2qn  +  eV  x  uT) 


or,  for  V  •  q*1  =  0  , 

5q n  =  -/i2V  x  (fj5jn  -  ew11) 

where 

wn  =  Vxf 


and 


|w"|  +  6 


r£(Rl+<r1i 

J _ 

N 


5.1.5 


where  the  sum  is  the  same  as  in  Eq.  5.1.3. 
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1.5.2  Results 

All  results  presented  are  in  2-D.  Also,  all  convection  cases  use  a  conservative 
second-order  centered  convection  term  together  with  a  centered  diffusion  term 
and  are  first-order-explicit  in  time.  In  all  the  result  plots,  axes  are  labeled 
with  grid  node  location.  Plots  of  amplitude  are  made  using  dense  contours 
with  white  grid  lines  superimposed. 

1.5. 2 .1  Passive  Scalar  Convection 

A  good  test  involves  convection  of  a  small  passive  scalar  pulse  in  solid-body 
rotating  flow(22>.  For  this  case  a  convection  step  was  followed  by  a  single 
confinement  step  (including  diffusion).  The  full  summation  of  Eq.  5.1.2  with 
Ci  =  1  (Eq.  5.1.3)  leads  to  a  spreading  of  the  pulses  in  the  direction  normal  to 
their  motion  (This  spreading  does  not  occur  for  vortices,  for  reasons  described 
in  Sec.  5.2.2).  This  scalar  spreading  seemed  to  be  due  to  the  inclusion  of 
“downwind”  values,  which  is  known  to  cause  instability.  By  setting  Ct  =  0  for 
downwind  points,  the  problem  is  almost  completely  cured,  and  only  a  very 
weak  spreading  remains. 

Results  are  shown  in  Fig.  7a  for  a  single  pulse  convecting  through  the  1st 
rotation  (2tt)  and  in  Fig.  7b  through  the  10th  rotation.  Contours  are  plotted 
from  initial  maximum  value  to  1/3  of  that  value.  The  final  image  in  Fig.  7b  was 
after  12,560  time  steps.  This  represented  a  travel  of  1,256  cells  or  about  433 
pulse  diameters.  The  open  circle  labels  the  starting  location.  The  time  step 
was  not  set  to  exactly  divide  the  circuit  so  we  would  not  expect  the  final  image 
to  exactly  coincide  with  the  initial  one.  However,  the  final  position  is  within 
plottable  accuracy  of  the  calculated  angle  and  radius.  The  final  radius  can  be 
seen  to  be  very  close  to  the  initial  one.  This  is  somewhat  surprising  since  the 
method  was  only  first  order  in  time  and  second  order  in  space.  We  have  not 
yet  attempted  to  reduce  the  small,  residual  cross-flow  spreading.  This  could 
easily  be  accomplished  with  a  small  amount  of  additional  negative  diffusion 
in  the  cross-stream  direction.  For  the  case  computed,  we  used  p  =  0.15  and 
e  =  0.75. 

Finally,  a  run  was  made  with  no  confinement  and  a  lower,  (minimal) 
diffusion  set  for  stability.  The  maximum  amplitude  decreased  to  below  the 
plotting  threshold  after  only  tt/8  radians  of  travel  (position  of  the  first  image 
in  Fig.  2a). 

1.5. 2. 2  Vorticity  Convection 

The  incompressible  fluid  dynamic  equations  with  Vorticity  Confinement, 
Eq.  3.0,  were  discretized  and  solved.  Each  time  step,  the  same  convection 
scheme  used  for  the  above  scalar  was  used.  A  standard  pressure  projection 
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method  was  used  to  enforce  incompressibility*23).  This  used  a  staggered  grid 
for  pressure  and  velocity  and  a  box  scheme  for  the  pressure  gradient.  The 
required  Poisson  equation  for  pressure  was  solved  using  a  direct  method 
“FishPack”.  Vorticity  Confinement  terms  -  Eqs  5.1.4  and  5.1.5,  were  used, 
with  n  =  0.05  and  e  =  1.5 fi.  Results  were  very  close  to  those  for  {i  —  0.2 
and  e  =  1.5 /x,  implying  that  the  method  is  not  very  sensitive  to  the  values  of 
these  parameters.  (Further  investigation  of  this  point  should  be  done).  The 
time  step,  based  on  grid  cell  size  and  circumferential  velocity  near  the  vortex 
cores,  was  0.4  and  the  grid  had  128  x  128  cells.  For  vortices,  the  full  summation 
(Eq.  5.1.3)  with  all  Ct  =  1  was  used,  without  excluding  downwind  points.  This 
exclusion,  useful  for  scalar  convection,  was  not  important  here  because  the 
large  circulating  velocities  near  the  vortex  cores  rapidly  convected  any  cross¬ 
stream  spread  vorticity  around  the  cores,  keeping  them  'close  to  axisymmetric. 
The  confinement  then  kept  the  cores  compact. 


First,  interacting  vortices  of  opposite  sign  were  computed  with  no  external 
free-stream- velocity,  convecting  only  under  their  own  induced  velocities.  To 
approximate  far-field  boundary  conditions,  each  time  step  the  total  amplitudes 
and  centroids  of  both  positive  and  negative  vorticity  were  computed.  The 
velocities  induced  by  corresponding  point  vortices  at  these  locations  were 
computed  and  imposed  on  the  boundaries.  Dirichlet  conditions  were  used  for 
pressure. 

Results  are  presented  in  Fig.  8  for  a  sequence  of  time  steps.  The  computed 
vortex  centroids  are  convecting  at  the  same  velocity,  to  plottable  accuracy,  as 
if  they  were  point  vortices.  The  self-induced  velocity  in  each  vortex,  although 
several  times  the  convecting  velocity  of  the  centroid,  automatically  cancels 
and  has  no  effect  because  the  method  conserves  momentum,  as  explained  in 
Sec.  5.1.1.  The  final  images  shown  result  from  1,800  time  steps,  yet  the  vortices 
are  just  as  compact  as  initially.  It  should  be  mentioned  that  if  the  vortices  are 
significantly  closer,  there  will  be  an  interaction  and  eventually  vorticity  will 
be  exchanged.  This  is  to  be  expected  since  there  are  “tails”  of  vorticity  that 
(rapidly)  decrease  with  radius  beyond  the  region  shown. 


Vortices  with  the  same  sign  were  computed  as  they  rotated  around  each 
other  under  their  induced  velocities.  The  vortices  were  initially  (and  finally) 
separated  by  14  cells.  First,  the  computation  was  done  with  no  Vorticity 
Confinement,  and  with  minimum  diffusion  required  for  convective  stability. 

Contours  are  displayed  in  Fig.  9  corresponding  to  vorticity  values  between 
initial  maximum  and  1  /4  of  that  value.  It  can  be  seen  that  the  vorticity  rapidly 
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diffuses.  The  computation  was  repeated  with  Vorticity  Confinement.  Values 
of  fi  and  e  were  those  used  above  for  opposite  sign  vortices.  Vorticity  contours 
between  initial  maximum  and  1/4  maximum  are  displayed  for  3  different  times 
in  the  1st  loop  around  each  other  in  Fig.  9a,  and  in  the  20th  loop  in  Fig.  9b. 
The  behavior  is  as  expected:  The  vortices  can  be  seen  to  be  essentially  the 
same  at  the  end  (after  7,200  time  steps)  as  initially. 

As  stated,  the  main  reason  for  using  an  Eulerian  method  as  opposed  to 
Lagrangian  is  to  allow  vortices  to  merge  and  otherwise  change  topology,  so 
that  complex,  realistic  vortical  configurations  can  be  treated.  This  has  been 
amply  demonstrated  by  the  original  Vorticity  Confinement  method  both  in 
2-D  and  3-D(4’5>7’9’10).  To  demonstrate  that  the  new  confinement  method 
does  not  inhibit  this  merging  we  repeat  the  above  computation,  but  with 
the  vortices  initially  5  cells  apart.  This  is  a  repeat  of  a  calculation  using  the 
original  Vorticity  Confinement  method(5>.  It  is  known  from  vorticity  contour  of 
“waterbag”  methods^24)  that  corotating  vortices  become  unstable  and  merge 
when  they  are  initially  closer  than  a  few  diameters.  It  can  be  seen  in  the 
vorticity  contour  plots  of  Fig.  10  that  this  merging  automatically  occurs,  as 
expected,  with  the  new  Vorticity  Confinement  method. 


1.6  Conclusion 

The  ability  of  the  original  Vorticity  Confinement  method  to  efficiently 
compute  blunt  body  turbulent  wakes  was  demonstrated.  This  involved  an 
LES-type  representation  but  with  the  small  scales  treated  as  inherently 
discrete  structures.  The  computations  involved  a  circular  and  square  cylinder 
in  3-D. 

A  new  Eulerian  Vorticity  Confinement  has  been  introduced.  Its  advantages 
over  the  original  Vorticity  Confinement  method  are  that  it  is  simpler  and  easier 
to  analyse  and  explicitly  conserves  total  momentum.  The  method  accurately 
convects  thin  features,  either  vortices  or  scalars,  over  very  long  distances, 
even  though  they  are  confined  to  only  2-3  grid  cells.  As  in  shock  and  other 
discontinuity  capturing  schemes,  the  details  of  the  internal  structure  of  the 
feature  are  not  solved  for,  but  implicitly  modeled  (the  features  are  treated 
as  weak  solutions).  Integral  quantities,  such  as  total  amplitude  and  centroid, 
however,  are  accurately  solved  for.  As  such,  the  method  essentially  represents 
the  features  as  nonlinear  discrete  solitary  waves  that  live  on  the  computational 
lattice.  It  is  argued  that  the  new  method  is  a  rational  generalization  of  1-D 
discontinuity  confinement  schemes  to  multiple  dimensions. 

Examples  are  presented  for  thin,  convecting  passive  scalars  and  vortices  in 
2-D  incompressible  flow. 
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(c)  ft  =  0.15,  e  =  0.  (No  Confinement)  (c)  (i  =  0.15,  e  =  0.  (No  Confinement) 

Figure  1  Mean  Streamwise  Velocity  Profiles  Figure  2  Streamwise  RMS  Fluctuations 
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(c)  \x  =  0.15,  e  =  0.0  (No  Confinement) 

Figure  3  Isosurfaces  of  Vorticity  Magnitude  (Isosurface  Level  =  0.15) 
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(b)  £  =  0.25 


Figure  4  Comparison  of  time-averaged  streamwise  velocity  along  a 
streamwise  line.  Symbols  denote  experimental  data. 


Figure  5  Comparison  of  time-averaged  velocity  profiles  at  x  =  1.  Symbols 
are  experimental  data. 
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(b)  €  =  0.25 


Figure  6  Comparison  of  root  mean  square  velocity  fluctuation  profiles  at 
x  =  1.  Symbols  are  experimental  data. 
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(a)  1st  loop  (b)  10th  loop 

Figure  7  Scalar  convection 
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Figure  8  Vorticity  contour  plots  of  self-induced  flow  by  two  point  vortices 
(15  cells  apart)  with  opposite  sign 
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Figure  9  Vorticity  contour  plots  of  self-induced  flow  by  two  vortices  (14  cells 
apart)  with  same  sign,  (a)  without  Vorticity  Confinement;  (b)with  Vorticity 
Confinement,  1st  loop;  (c)with  Vorticity  Confinement,  after  20th  loops. 
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Abstract 


Over  the  last  few  years,  a  new 
flow  computational  methodology, 
Vorticity  Confinement,  has  been  shown 
to  be  very  effective  in  treating 
concentrated  vortical  regions.  These 
include  thin  vortex  filaments  which  can 
be  numerically  convected  over  arbitrary 
distances  on  coarse  Eulerian  grids  while 
requiring  only  ~2  grid  cells  across  their 
cross  section.  They  also  include 
boundary  layers  on  surfaces  “immersed” 
in  non-conforming  uniform  Cartesian 
grids,  with  no  requirement  for  grid 
refinement  or  complex  logic  near  the 
surface. 

In  this  paper  we  use  Vorticity 
Confinement  to  treat  flow  over  blunt 
bodies,  including  attached  and 
separating  boundary  layers,  and  resulting 
turbulent  wakes.  The  same  basic  idea  is 
applied  to  all  of  these  features:  At  the 
smallest  scales  (~2  cells)  the  vortical 
structures  are  captured  and  treated, 
effectively,  as  solitary  waves  that  are 
solutions  of  nonlinear  discrete  equations 
on  the  grid.  These  do  not  represent 
discretizations  of  partial  differential 
equations  (pde’s)  but  nonlinear  models 


of  the  structures,  directly  on  the  grid. 
They  allow  the  boundary  layer  to  be 
effectively  “captured”.  In  the  turbulent 
wake,  where  there  are  many  scales,  they 
represent  an  effective  small  scale  energy 
sink.  However,  they  do  not  have  the 
unphysical  spreading  due  to  numerical 
diffusion  at  these  scales,  which  is 
present  in  conventional  computational 
methods. 

The  basic  modeling  idea  is 
similar  to  that  used  in  shock  capturing, 
where  intrinsically  discrete  equations  are 
satisfied  in  thin,  modeled  regions.  It  is 
argued  that,  for  realistic  high  Reynolds 
number  flows,  this  direct,  grid-based 
modeling  approach  is  much  more 
effective  than  first  formulating  model 
pde’s  for  the  small  scale,  turbulent 
vortical  regions  and  then  discretizing 
them. 

Results  are  presented  for  3-D 
flows  over  round  and  square  cylinders 
and  a  realistic  helicopter  landing  ship. 
Comparisons  with  experimental  data  are 
given. 
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Finally,  a  new  simpler 
formulation  of  Vorticity  Confinement  is 
given  together  with  a  related  formulation 
for  confinement  of  passive  scalar  fields. 


1.  Introduction 

Many  real-world  flows  are 
characterized  by  regions  of  concentrated 
vorticity.  These  are  typically  turbulent  at 
realistic  Reynolds  numbers  and  can 
convect  over  long  distances,  either  as 
thin  vortex  filaments  or  blunt  body 
wakes  containing  small-scale  vortical 
structures.  Also,  boundary  layers  are 
typically  thin,  turbulent  vortical 
structures  near  solid  surfaces.  Even 
though  the  Navier-Stokes  equations 
apply  to  these  flows,  it  is  not  feasible  to 
solve  them  on  foreseeable  computers 
due  to  the  ever-present  small  scales.  As  a 
result,  conventional  CFD  methods 
involve  first  formulating  the  partial 
differential  equations  (pde’s)  that  model 
these  turbulent  regions.  These  pde’s  are 
then  discretized  in  the  turbulent  regions 
and  then  solved.  This  typically  requires 
large  computer  resources  and  difficult 
grid  generation.  The  problem  is  that 
resolving  even  the  model  pde’s  requires 
very  fine  computational  grids,  which 
must  conform  to  the  surface  for  the 
boundary  layer.  Further,  these  methods 
typically  dissipate  thin  filaments  and 
vortical  structures  in  the  wake  as  they 
convect.  This  is  true  even  if  on  the  order 
of  10  grid  points  are  devoted  to  the  cross 
section  of  each  structure  and  complex, 
high  order  discretizations  are  used.  The 
problem  is  intrinsic  to  the  discretization 
of  the  convective  terms  and  made  worse 
by  the  use  of  dissipative  terms  to  model 
the  turbulent  effects.  Finally,  even  if  the 
model  pde’s  are  resolved,  the  result  is, 


of  course,  still  a  model  approximation  of 
the  turbulent  flow. 

In  this  paper,  we  describe  an 
entirely  new  CFD  methodology  where 
very  different  numerical  models  are  used 
for  vortical  regions.  Instead  of  first 
hypothesizing  a  turbulence  model  based 
on  pde’s  and  then  trying  to  accurately 
discretize  and  resolve  them,  we  model 
the  internal  structure  of  the  vortical 
regions  directly  on  the  grid  using 
generalized  nonlinear  difference 
equations,  rather  than  using  finite 
difference  discretizations  of  partial 
differential  equations  that  attempt  to 
approximately  resolve  the  model 
equations  (and  the  vortical  structures). 
This  approach  allows  treatment  of 
vortical  structures  as  objects  spread  over 
only  1-2  grid  cells  on  coarse,  essentially 
uniform  Cartesian  computational  grids. 

These  ideas  are  implemented  in 
the  methodology  termed  Vorticity 
Confinement.  Although  new  for  vortical 
regions,  these  ideas  are  old  for  the 
treatment  of  shocks,  starting  with 
VonNeumann  &  Richtmyer  [1],  Lax  [2] 
and  others.  There,  as  is  well  known, 
shocks  are  treated  as  thin  regions  spread 
over  a  few  grid  cells  that  obey  discrete, 
grid-based  nonlinear  model  equations, 
that  conserve  certain  quantities.  This 
approach  has  been  shown  to  be  more 
effective  than  trying  to  discretize  and 
solve  the  applicable  Navier-Stokes  pde’s 
in  those  thin  regions.  One  difference, 
however,  is  that  with  shocks,  unlike 
vortical  regions,  characteristics  slope 
inward  toward  the  shock.  As  a  result,  the 
modeling  is  simpler  than  that  with 
vortical  structures.  One  of  the  first 
confinement-type  schemes  for  contact 
discontinuities,  where,  unlike  shocks, 
physical  characteristics  do  not  help,  was 
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developed  by  Harten  [3],  but  was 
specialized  to  one-dimensional 
compressible  flow. 

Many  results  have  been  obtained 
in  recent  years  using  Vorticity 
Confinement  [4-22].  Lohner  has  a  very 
good  short  review  of  this  work  [18]. 
Initial  results  were  for  isolated, 
convecting  vortex  filaments  [14-16]. 
Then,  vortices  convecting  past  airfoils 
and  wings  (blade  -  vortex  interactions) 
were  treated  [13].  In  this  early  study, 
unlike  in  our  current  studies,  near  the 
surface  a  surface-fitted  grid  was  used  for 
the  wing  with  surface  grid  refinement  to 
resolve  the  actual  Navier-Stokes 
equations,  since  only  a  low  Reynolds 
number,  laminar,  case  was  treated.  To 
accommodate  this  grid  refinement  with 
Vorticity  Confinement,  the  parameter 
specifying  the  strength  of  the  Vorticity 
Confinement  term  (s  )  was  made  to  be 
proportional  to  grid  size  so  that  it 
automatically  vanished  in  the  fine-grid 
boundary  layer  region,  but  was  able  to 
confine  the  convecting  vortex  in  the 
external,  coarse-grid  region.  Our  current 
studies,  like  those  described  in  this 
paper,  involve  surfaces  “immersed”  in 
uniform,  non-conforming  grids  with  no 
grid  refinement  and  use  a  constant  value 
for  £  . 

Recently,  Vorticity  Confinement 
has  been  used  together  with  unstructured 
grids  [18,  19].  When  using  these  grids, 
which  have  rapidly  changing  cell  sizes, 
care  must  be  taken  not  only  that  £  varies 
properly  with  cell  size,  but  also  that  the 
confinement  correction  does  not  extend 
beyond  the  vortex  core  due  to  numerical 
artifacts  of  the  implementation.  This 
property  is  true  in  the  continuum  limit, 
as  shown  in  the  description  of  Sec.2,  and 
should  be  preserved  in  the  discretization. 
If  the  correction  does  extend  beyond  the 


vortex,  then  it  could  erroneously  affect 
surface  pressure  if  a  vortex  is  passing 
near  a  surface.  This  could  be  important, 
for  example,  for  delta  wings  [18]  and 
similar  cases,  where  vortices  convect 
near  surfaces.  In  fact,  for  delta  wings,  it 
is  well  known  that  there  is  a  feeding 
sheet  from  the  leading  edge  causing  the 
vortex  to  grow  in  strength  as  it  convects 
and  causing  the  characteristics  to  point 
towards  it.  In  such  cases,  for  a 
reasonable  grid,  confinement  is  not 
really  needed  (until  the  vortex  convects 
past  the  trailing  edge).  If  confinement  is 
used  correctly,  however,  it  should  not 
change  the  nearby  pressure  on  the 
surface  even  in  these  cases,  for  high 
Reynolds  number  flow.  Finally,  in  low 
Reynolds  number  viscous  flow  cases 
where  laminar  flow  occurs,  the  proper 
Navier-Stokes  equations  should,  of 
course,  be  used,  since  the  vortical  scales 
are  then  not  small  and  there  is  no  need  to 
discretize  model  equations  for  small 
scales. 

In  this  paper,  Vorticity 
Confinement  is  applied  to  a  series  of 
blunt  bodies,  a  circular  and  a  square 
cylinder,  and  a  realistic  helicopter 
landing  ship.  In  the  cylinder  cases,  for 
which  unsteady  experimental  results  are 
available,  comparison  is  made.  The  ship 
case  demonstrates  the  ability  of  Vorticity 
Confinement  to  preserve  thin, 
concentrated  vortical  structures  over 
long  distances.  Also  velocity  and  vortex 
trajectory  computations  are  shown  to 
agree  well  with  experimental  results. 
All  cases  were  run  on  grids  much 
coarser  than  those  required  by 
conventional  CFD  methods  and  did  not 
require  body  conforming  grid  generation 
or  refinement  near  the  surfaces. 
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2.  Vorticitv  Confinement 

Vorticity  Confinement  is  a 
method  to  preserve  convecting, 
concentrated  vorticity  on  a  coarse  grid; 
the  vortical  structures  can  represent 
boundary  layers  over  solid  surfaces,  the 
smaller  scales  (of  the  order  of  a  grid  cell) 
in  turbulent  wakes,  or  separated,  thin 
vortex  filaments  that  convect  over 
arbitrarily  long  distances.  Vorticity 
Confinement  can  be  implemented  in  a 
pre-existing  flow  solver,  for  both 
incompressible  and  compressible  flow, 
by  adding  a  term  to  the  discretized 
momentum  conservation  equations  [12]. 

For  general  unsteady 
incompressible  flows,  the  governing 
equations  with  the  Vorticity 
Confinement  term  are  discretizations  of 
the  continuity  equation  and  the 
momentum  equations,  with  an  added 
term: 

V-<7  =  0  2.0.1 

d,q  =  -(q-V)q-— Vp  +  [pV2^-£5] 

P 

2.0.2 

where  q  is  the  velocity  vector,  p  is  the 
pressure,  p  is  the  density,  and  p  is  a 
diffusion  coefficient  that  includes 
numerical  effects  (we  assume  physical 
diffusion  is  much  smaller).  For  the  last 
term,  es ,  e  is  a  numerical  coefficient 
that,  together  with  p ,  controls  the  size 
and  time  scales  of  the  convecting 
vortical  regions  or  vortical  boundary 
layers.  For  this  reason,  we  refer  to  the 
two  terms  in  the  brackets  as 
“confinement  terms”. 


The  first  confinement  term,  or 
diffusion  term  is  usually  specified 
explicitly,  but  can  be  implicitly  present 
in  the  solver,  such  as  in  a  lower  order 
upwind  scheme.  There  are  also  many 
possible  forms  for  the  second 
confinement  term.  First,  the  original  one 
used  in  this  and  earlier  studies  will  be 
described.  Then,  a  new  more  elegant, 
simpler  form  will  be  described,  together 
with  a  simple  demonstration. 


2.1  Original  Form 

S  =  «X( 0 

2.1.1 

where 

n  =  Vq  /|Vt|| 

2.1.2 

the  vorticity  vector  is  given  by 

(5  =  Vxq 

2.1.3 

and 

rl=N 

2.1.4 

In  general,  for  boundary  layers 
and  convecting  vortex  filaments, 
computed  flow  fields  external  to  the 
vortical  regions  are  not  sensitive  to  the 
parameters  £  and  p  over  a  wide  range 
of  values.  For  example,  the  flow  outside 
an  axisymmetric  2-D  vortex  core  is 
independent  of  the  vortical  distribution, 
and  hence  does  not  depend  on  e  and  p 
as  long  as  the  core  is  thin.  Hence,  the 
issues  involved  in  setting  them  are 
similar  to  those  involved  in  setting 
numerical  parameters  in  other  standard 
computational  fluid  dynamics  schemes, 
such  as  artificial  dissipation  in  many 
conventional  compressible  solvers. 
There,  the  flow  outside  1-D  captured 
shock  regions  does  not  depend  on  the 
exact  internal  structure,  as  long  as  it  is 
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thin.  For  wake  flows,  e  and  p  can  be 
used  to  approximately  simulate  finite 
Reynolds  number  effects,  since  they 
control  the  intensity  of  the  small  vortical 
scales. 

An  important  feature  of  the 
Vorticity  Confinement  method  is  that  the 
Confinement  terms  are  non-zero  only  in 
the  vortical  regions,  since  both  the 
diffusion  term  and  the  anti-diffusion 
term  vanish  outside  those  regions  (care 
has  to  be  taken  in  the  numerical 
implementation  to  preserve  this  feature). 

Another  important  feature 
concerns  the  total  change  induced  by  the 
confinement  correction  in  mass,  vorticity 
and  momentum,  integrated  over  a  cross 
section  of  a  convecting  vortex.  It  can  be 
shown  [14,16]  that  mass  is  conserved 
because  of  the  pressure  projection  step  in 
the  solver  and  vorticity  is  explicitly 
conserved  because  of  the  vanishing  of 
the  correction  outside  the  vortical 
regions.  Momentum  is  almost  exactly 
conserved.  A  small  extension  of  the 
original  method  [9],  explicitly  conserves 
the  momentum.  The  new  formulation, 
described  below  in  Sec.2.2  also 
explicitly  conserves  momentum. 

It  should  be  mentioned  that  the 
above  pde  form  for  the  confinement 
cannot  be  said  to  be  resolved  since 
vorticies  are  captured  in  only  a  few  grid 
cells.  Accordingly,  we  can  only  say  that 
the  pde  form  is  only  used  to  motivate 
the  final  discrete  form  and  that  a  large 
body  of  numerical  evidence  has  been 
given  that  the  terms  result  in  confined 
vortices. 


Formulation 


Because  of  the  above  issue 
concerning  the  pde  form  and  that 
additional  corrections  must  be  added  to 
make  it  explicitly  conserve  momentum, 
a  new,  simpler  formulation  that  does  not 
have  these  problems  has  been  developed. 
This  new,  intrinsically  discrete 
formulation  is  presented  in  this  section. 

A  more  detailed  description  is 
being  presented  in  [23].  First,  a 
formulation  for  scalar  confinement  is 
given.  Then,  a  velocity-based  Vorticity 
Confinement  correction  is  given  that 
reduces  to  the  scalar  confinement  in 
terms  of  vorticity  when  the  curl  is  taken. 

The  scalar  formulation  presented 
here  is  related  to  that  presented  in  [12]  in 
1-D.  In  this  section  no  convection  is 
used,  only  the  two  confinement  terms,  so 
that  the  behavior  can  be  seen  more 
clearly.  Excellent  results  are  found  with 
convection  and  will  be  shown  in  [23]. 

We  start  with  an  iteration  for  a 
scalar  4» : 

<j>n+1  =<|)"  +  ^VV-£V20"  2.2.1 
where 

XcA'+s)-T' 

<D"  =  - — = -  2.2.2 

ic, 

/ 

where  the  sum  is  over  a  set  of  grid  nodes 
near  and  including  the  node  where  O  is 
computed,  and  8  is  a  small  positive 
constant  (~  10“8  )  to  prevent  problems 
due  to  finite  precision.  The  coefficients, 
C,,  can  be  varied  but  good  results  are 

obtained  by  simply  setting  them  to  1. 

For  example,  in  2-D,  one 
possibility  is 
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where  N  neighboring  terms  are  taken. 

Here,  we  assume  <|) "  >  0 . 
Negative  values  can  also  be 
accommodated  with  a  small  extension. 
Both  (x  and  £  are  positive. 

An  important  feature  is  that  all 
terms  are  homogeneous  of  degree  1  in 
Eq.  2.2.1.  This  is  important  because  the 
confinement  not  depend  on  the  scale  of 
the  quantity  being  confined.  Another 
important  feature  is  the  nonlinearity  of 
different  order  in  the  derivatives.  It  is 
easy  to  show  that  a  linear  combination  of 
terms  cannot  lead  to  a  stable 
confinement  for  any  finite  range  of 
coefficients. 

For  smooth  <|)  fields,  the  last 
term  represents  a  negative  diffusion.  If 
p  <  e  ,  where  N  is  the  number  of  terms 
in  Eq.  2.2.2,  the  total  diffusion  is 
negative.  However,  the  iteration  Eq. 
2.2.1  is  still  stable  and  converges  for 
values  of  £  up  to  about  5  p ,  resulting  in 
an  effective  negative  diffusion 
coefficient  in  the  long  wavelength  limit. 
Also,  the  discrete,  converged  solution 
can  be  given  exactly  in  terms  of  sech 
functions. 

The  next  step  involves  letting  <j> 
be  the  magnitude  of  vorticity  and 
deriving  an  equation  for  the 
corresponding  velocity  correction  that 
leads  to  Eq.  2.2.1  when  the  curl  is  taken 
(exactly  in  2-D,  or  for  a  straight  vortex 


in  3-D).  This  will  result  in  a  new 
formulation  of  Vorticity  Confinement. 

We  simply  define 

qn+l  =  q"  +  pV2tf"  +eVxwn  2.2.4 

or,  for  V  •  q  =  0  , 

qn+x  =qn  -Vx(p(o"-£w")  2.2.5 

where 

<5  n=Vxq"  2.2.6 


and 


where  the  sum  is  the  same  as  in  Eq. 

2.2.2. 

Results  are  presented  in  Fig.l 
after  0,  8,  and  100  iterations,  for 
vorticity  and  velocity  for  a  vortex  in  2- 
D.  Values  used  were  p  =  .2 ,  £  =  0  and 
a=5p.  In  this  figure,  the  vorticity 
contour  levels  extend  from  about  !4  of 
the  maximum  initial  value  to  the 
maximum  so  that  a  measure  of  the  size 
of  the  confined  region  can  be 
determined.  We  are  currently 
investigating  the  new  method  for 
implementation  in  our  codes.  In 
addition,  the  use  of  the  “Scalar 
Confinement”  version,  Eq.  2.2.1,  is 
being  used  for  flows  where  thin  streams 
of  passive  scalars,  such  as  contaminants, 
must  be  convected  over  long  distances. 
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3.  Results 


a)  initial  condition 


b)time  steps  =  8  £  =  0 


d)  time  steps  =  8 

£  =  5  |_l 


e)  time  steps  =100 
£  =  5(J. 

Figure.  1  Vorticity  contours  and 
vector  fields  of  velocity 


Three  test  cases  are  presented  in 
this  section.  First,  the  effects  of 
Vorticity  Confinement  on  the  wake 
structure  of  a  3-D  circular  and  square 
cylinder  are  presented.  The  third 
example  shows  the  application  of 
Vorticity  Confinement  to  a  realistic  flow 
simulation  for  a  ship  configuration. 

3.1  3-D  Cylinder 
3.1a  3-D  Circular  Cylinder 

Flow  over  a  3-D  circular  cylinder 
was  calculated  to  assess  the  ability  of 
Vorticity  Confinement  to  accurately 
model  the  wake  flow  behind  a  blunt 
body.  A  long  cylinder  was  “immersed” 
in  a  uniform  141x101x61  Cartesian 
grid  in  the  streamwise,  normal  and 
spanwise  directions,  respectively. 
Periodic  conditions  were  imposed  at  the 
lateral  boundaries.  The  diameter  of  the 
cylinder  was  15  grid  cells.  The  origin  of 
the  coordinate  system  used  is  located  in 
the  center  of  the  cylinder,  and  all 
distances  are  non-dimensionalized  by  the 
diameter,  shown  in  Fig.  2.  Results  of 
the  computations  were  compared  to  the 
experimental  results  of  Lourenco  and 
Shih  at  a  Reynolds  number  of  about 
3,900,  described  in  Ref.  [24].  The 
computational  results  were  all  averaged 
over  the  spanwise  direction. 

The  diffusion  coefficient  p  was 
held  constant  for  this  study.  The 
confinement  coefficient,  s  ,  was  adjusted 
to  impose  different  levels  of 
confinement.  This  resulted  in  different 
levels  of  the  intensity  of  the  small 
vortical  scales  in  the  wake,  and 
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approximately  simulated  different 
Reynolds  numbers. 


Figure  2.  Measurement  Positions  For 
Circular  Cylinder 

Figure  3  depicts  the  mean 
streamwise  velocities  resulting  from 
three  values  of  the  confinement 
coefficient.  Figures  a)  and  b)  depict  the 
result  of  two  levels  of  confinement, 
where  the  flow  is  measured  along  lines 
which  are  normal  to  the  cylinder  axis 
and  the  mean  stream,  at  three  different 
locations  in  the  wake  of  the  cylinder 
shown  in  Fig.  2.  For  both  values  of  e 
the  agreement  can  be  seen  to  be  very 
good,  indicating  that  the  effect  of  the 
confinement  parameter  is  small  over  a 
range  of  values.  Results  from  a  case 
without  confinement  are  depicted  in  Fig. 
3  c.  Without  confinement,  the  flow  field 
is  dominated  by  diffusive  effects  that  are 
not  counterbalanced  by  the  anti-diffusive 
confinement  term  and  approximates  a 
steady,  low  Reynolds  number  flow. 


c)  p  =  0.15,  e  =  0.0  (No  Confinement) 

Figure  3.  Mean  Streamwise  Velocity 
Profiles.  Symbols  are  experimental  Data. 


The  ability  of  Vorticity 
Confinement  to  model  a  turbulent  wake 
was  assessed  by  computing  the  rms 
streamwise  velocity  fluctuations  in  the 
wake  region.  Comparisons  of  these 
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fluctuations  with  the  experimental  data 
were  made  along  the  same  lines  in  the 
wake  as  for  the  mean  velocity,  for  the 
same  three  values  of  e  .  Fig.  4a  (s  = 
0.25)  shows  very  good  agreement.  The 
effect  of  increasing  confinement  to  0.5  is 
depicted  in  Fig.  4b.  In  general,  the 
effect  of  increased  confinement  is  to  thin 
the  shear  layer  comprising  the  wake 
boundary,  and  to  increase  the  fluctuation 
of  the  time-dependent  flow  from  the 
mean  flow.  Fig.  4c  depicts  results 
without  confinement.  Without 

confinement,  the  flow  can  be  seen  to  be 
steady  and  exhibits  none  of  the 
fluctuations  that  would  occur  when 
confinement  is  used. 

Figure  5  depicts  isosurfaces  of 
vorticity  magnitude  for  the  same  three 
levels  of  confinement.  The  use  of 
confinement  results  in  chaotic  flow 
patterns,  as  would  be  expected  in  three- 
dimensional  turbulent  flows.  This  does 
not  occur  in  two-dimensional 
simulations,  indicating  that  this  chaotic 
behavior  is  not  due  to  numerical 
instability  created  by  the  confinement. 
Increasing  confinement  increases  the 
chaotic  nature  of  the  flow  in  3-D  and 
reduces  the  characteristic  size  of  the 
vortical  structures,  analogous  to  what 
would  be  expected  with  an  increase  in 
Reynolds  number.  Clearly,  the  use  of 
confinement  allows  small-scale  time- 
dependent  wake  structures  to  be 
generated  on  extremely  coarse  grids  with 
the  small-scale  structure  captured  over 
only  1~2  grid  cells.  Also,  these  small- 
scale  structures  serve  as  a  viscous  sink 
for  turbulent  energy,  as  in  physical 
turbulence. 


a)  p  =  0.15,  8  =  0.25 
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c)  p  =  0.15,  e  =  0.0  (no  confinement) 

Figure  4.  Streamwise  Reynolds 
Stresses.  Symbols  Denote  Experimental 
Data 
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a)  p  =  0.15,  e  =  0.25,  Isosurface  Level 
=  ±0.15 


b)  p  =  0.15,  e  =  0.5,  Isosurface  Level  = 
±0.15 


c)  p  =0.15,  e  =  0.,  Isosurface  Level  = 
±0.15 

Figure  5.  Isosurfaces  of  Vorticity 
Magnitude 


3.1b  3-D  Square  Cylinder 

Flow  over  a  square  cylinder  was 
also  calculated.  As  in  the  circular  case, 
the  cylinder  was  “immersed”  in  a 
uniform  141x101x61  Cartesian  grid 
and  periodic  conditions  imposed  at  the 
lateral  boundaries.  The  diameter  (length 
of  each  side)  of  the  cylinder  was  also  15 
grid  cells.  The  same  coordinate  system 
was  used  as  for  the  circular  cylinder,  as 
shown  in  Figure  6.  Results  of  the 
computations  were  compared  to  the 
experimental  results  of  Lyn  et  al.  [25]  at 
a  Reynolds  number  of  about  21,400.  The 
computational  results  were  all  averaged 
over  the  spanwise  direction. 


U 


XtD  *  1 

Figure  6:  Measurement  Positions  For 
Square  Cylinder 

As  in  the  circular  cylinder  case, 
the  diffusion  coefficient  p  was  also  held 
constant  at  0.15.  The  confinement 
coefficient,  e  ,  was  adjusted  to  impose 
different  levels  of  confinement  so  as  to 
approximate  the  effects  of  different 
Reynolds  numbers. 

Figure  7  depicts  the  comparison 
with  experimental  data  of  the  time- 
averaged  streamwise  velocity  along  a 
streamwise  line  extending  downstream 
from  the  middle  of  the  leeward  face  of 
the  cylinder.  Results  of  two  values  of  the 
confinement  coefficient  are  plotted. 
Figure  8  shows  the  time-averaged 
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velocity  along  a  line  normal  to  the 
cylinder  axis  and  the  mean  stream  at  x  = 
1.  Symbols  represent  the  experimental 
data.  Our  numerical  results  agree  well 
with  the  experimental  data. 


(a)  e  =0.35 


(b)  e  =  0.25 

Figure  7.  Comparison  of  time-averaged 
streamwise  velocity  along  a  streamwise 
line.  Symbols  denote  experimental  data. 

Comparisons  of  the  computed 
turbulence  level  with  the  experimental 
results  also  show  reasonably  good 
agreement  in  Fig.  9a.  The  effect  of 
decreasing  confinement  is  shown  in  Fig. 
9b.  Comparing  this  case  with  the 
previous  3-D  circular  cylinder  case,  it  is 
easy  to  see  that  a  higher  value  of  e  is 
required  for  better  results.  This  is 
apparently  because  the  Reynolds  number 
has  increased  from  about  3,900  in  the  3- 
D  circular  cylinder  case  to  roughly 
21,400  in  the  present  3-D  square 
cylinder  case. 


(a)  e  =0.35 


1J5, 


2  if  1  0*  0 

y 

(b)  e  =0.25 

Figure  8:  Comparison  of  time-averaged 
velocity  profiles  at  x  =  1.  Symbols  are 
experimental  data. 
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3.1c  Synopsis  of  Cylinder  Study 


(a)  s  =  0.35 


(b)  e  =0.25 

Figure  9:  Comparison  of  root  mean 
square  velocity  fluctuation  profiles  at 
x  =  1 .  Symbols  are  experimental  data. 


In  subsequent  studies,  the 
correlation  between  e  and  Reynolds 
number  will  be  studied  by  comparisons 
between  computation  and  experiment  for 
other  cases  and  a  useful  calibration  of  e 
determined.  Also,  it  should  be 
emphasized  that  our  agreement  with 
experiment  is  closer  than  many  much 
finer  grid  studies  with  much  more 
complex  LES  pde  models  [26],  each 
requiring  a  number  of  empirical 
coefficients.  Just  one  coefficient  is 
adjusted  here,  on  a  simple,  coarse  grid. 

3.2  Ship  Configuration 

The  present  method  of 
developing  the  operating  limits  for 
helicopters  landing  or  taking  off  from 
ships  is  accomplished  largely  through 
flight  tests,  which  are  time-consuming, 
costly,  and  potentially  dangerous.  The 
ability  to  develop  the  operational  limits 
using  computational  methods  as  an 
adjunct  to  present  methods  has  the 
potential  of  significantly  reducing  cost, 
time,  and  risk.  Figure  10  depicts  the  ship 
configuration  used  in  this  study,  which 
has  undergone  wind  tunnel  testing  by 
Caradonna  [27]. 

The  ship  configuration  is  80  units 
from  bow  to  stem.  The  front  of  the  bow 
is  located  at  X=0;  the  Y=0  plane  is 
perpendicular  to  the  deck  and  is 
coincident  with  the  deck  midplane.  The 
deck  is  13  units  wide,  extending  from 
Y—6.5  to  6.5.  The  deck  surface  is 
located  at  Z=0,  with  positive  Z  upwards. 
The  aerodynamic  domain  was  modeled 
with  two  grids:  a  401  x  121  x  74  grid 
extended  from  X=-50  to  150,  Y=+/-  30, 
and  Z=-6.5  to  30.  Far-field  boundary 
conditions  were  imposed  on  the  surfaces 


12 


of  this  grid.  A  second  grid  was  used  to 
more  highly  resolve  the  forward  deck 
region.  The  second  grid  (501  x  92  x  51) 
extended  from  X=-10  to  50,  Y=-8  to 
2.92,  and  Z=-2.5  to  3.5.  Boundary 
conditions  were  imposed  on  the  inner 
grid  via  interpolation.  In  addition,  the 
outer  grid  obtained  flow  field 
information  from  the  inner  grid  via 
interpolation  at  interior  boundaries.  As  a 
result,  two-way  flow  field 
communication  was  effected  to  obtain  a 
globally  consistent  flow  field.  Velocities 
were  measured  on  a  plane  (X=23.63)  for 
a  number  of  Z  locations  and  were 
compared  with  computational  results. 


Figure  10.  Helicopter  Landing  Ship 

Figure  1 1  depicts  vorticity 
isosurfaces  over  the  ship  deck  for  a  wind 
aspect  angle  of  20°.  Ship  flow  fields  are 
often  characterized  by  the  development 
of  a  strong  vortex  at  the  windward  edge 
of  the  deck  that  subsequently  convects 
across  the  ship  deck.  The  effect  of  field 
confinement  is  illustrated  in  a 
comparison  between  the  isosurfaces  of 
Figs.  1  la  and  lib.  Without  confinement 
(Fig.  11a),  the  vortex  is  greatly 
dissipated;  with  confinement  (Fig.  lib), 
the  deck  vortex  persists  indefinitely  over 
the  deck  surface.  The  view  in  Fig.  11  is 


looking  down  from  above  the  deck  with 
the  bow  at  the  bottom  of  the  figure. 


a)  No  Confinement  b)  With  Confinement 

Figure  11.  Isosurface  of  Vorticity  on 
Ship  Deck 

The  computed  vortex  location  is 
very  close  to  the  experimental  location. 
These  vortical  structures  are  major 
features  of  the  ship  flow  field,  and  would 
be  difficult  to  resolve  and  maintain  with 
conventional  CFD  methods  without  a 
much  finer  mesh.  By  contrast,  the 
vortical  structures  are  resolved  and 
persist  on  a  relatively  coarse  mesh  with 
Vorticity  Confinement.  This  is  in 
accordance  with  physical  expectations, 
in  that  the  strength  and  structure  of  the 
vortices  are  maintained  without 
significant  dissipation.  Quantitative 
results  are  presented  in  Fig.  12,  which 
depicts  velocity  one  unit  above  the  deck 
at  a  constant  height  horizontal  line  that 
goes  across  the  deck  through  the  vortex 
center.  The  line  is  depicted  (from  above) 
in  Fig.  lib  as  the  “experimental  plane”. 
Additional  information  on  the 
computation  is  given  in  Ref.  [4]. 
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a)  Horizontal  Velocity 


b)  Vertical  Velocity 

Figure  12.  Comparison  of  Experimental 
and  Computed  Velocity  Over  Deck, 
Z=1.0 

4.  Conclusions 

Vorticity  Confinement  has  been 
shown  to  provide  an  efficient,  promising 
means  of  computing  the  turbulent  wake 
behind  a  round  and  a  square  cylinder,  on 
a  veiy  coarse  mesh.  Results  were 
presented  for  the  mean  flow  and  the  rms 
fluctuations,  and  the  effects  of  different 
levels  of  confinement  (and  no 
confinement)  were  demonstrated.  It 
appears  that  the  ability  to  resolve  scales 


down  to  1~2  grid  cells  and  absorb 
energy  there  suffices  to  quantitatively 
solve  for  the  dynamics  of  the  large  scale 
eddies  in  the  wake.  Also,  the  method 
was  shown  to  allow  the  emdedding  of 
the  surface  in  a  uniform,  coarse 
Cartesian  grid  without  body-fitting, 
refinement,  or  complex  logic  near  the 
surface. 

It  can  be  seen  that  the  Vorticity 
Confinement  results  for  the  cylinder 
wakes,  when  compared  to  experiment, 
are  comparable  to  conventional  LES 
results.  The  difference  is  that  the 
confinement  solutions  employ  immersed 
boundaries  with  simple,  coarse  Cartesian 
grids,  with  only  16  grid  cells  across  the 
cylinder  diameter,  and  no  complex  LES 
turbulence  pde  models.  As  a  result,  no 
grid  generation  was  required  and  the 
computations  were  very  simple  to  set  up 
and  fast  to  run.  Based  on  these  results,  it 
appears  that  Vorticity  Confinement 
should  be  very  useful  for  rapid 
engineering  solutions  for  complex  flows. 
Additional  validation  is  being  done  for 
surface  pressures  and  forces  in  these 
cases  (beyond  that  already  done)  [5,  9]. 
Further,  the  functional  dependence  of 
our  single  constant,  s  ,  on  the  Reynolds 
number  is  being  calibrated  by  computing 
more  cases. 

Vorticity  Confinement  was  also 
applied  to  a  helicopter  landing  ship.  In 
this  case,  the  use  of  confinement  was 
shown  to  be  necessary,  for  the  grid  used, 
for  the  development  of  a  deck  vortex. 
This  vortex  is  a  major  feature  of  ship 
flow  fields  at  general  wind  aspect 
angles.  The  ship  solution  is  of  particular 
interest  because,  unlike  the  cylinder,  the 
main  vortical  structures  are  not  shed,  but 
remain  attached  to  the  configuration. 
Although  a  very  coarse  mesh  was  able  to 
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reproduce  small-scale  turbulent  structure 
quite  well  with  the  cylinder,  the  ship 
required  a  more  dense  mesh  near  the 
corner  of  the  bow  to  correctly  predict  the 
vortex  generation.  Similar  results  have 
been  found  in  other  cases.  This  is  due  to 
the  interaction  of  the  confinement  terms 
with  boundary  conditions  imposed  by 
solid  surfaces.  This  interaction  is  a 
subject  of  current  research,  as  many 
structures  (e.g.,  buildings)  may  exhibit 
bound  vortices,  and  the  ability  to  model 
these  structures  with  coarse  grids  is  of 
great  importance. 

Finally,  a  new,  more  elegant, 
fully  conservative  Vorticity 

Confinement  formulation  was  presented 
with  preliminary  results  with  diffusion 
but  without  convection.  This  will  be 
used  in  the  future. 
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